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16  SUPPLEMENTARY  NOTATION 


COSATI  COOES 


IS  SUBJECT  TERMS  (Continue  on  reverte  if  neceuary  and  i^nUfyby 

Attenuation,  Magnitude,  Seismic  Yield,  Q,  Lg  Waves, 
Surface  Waves,  Spectra,  Nuclear  Explosions 


’’  "“Tre'VbTecT/v^r^  explosions  are  to  obtain  ^Ji^ration 

curves  or  equations  of  mb(Lg)  versus  the  logarithm  of  the  explosion  yield  and  to  obtain 
estimates  of  the  mb(P)  bias  between  test  sites,  such  as  NTS  and  Shagan  River, 
nf  mKfLel  values  The  first  objective  will  provide  a  means  of  estimating  directly  t 
e^Iolion  TieirfrorLfsured  L^g-wave  amplitudes.  The  second  objective  wHl 
m^P)  bias  information  that  is  needed  to  estimate  yields  from  mb(P)  values  at  places 

other  than  NTS. 

Technical  problems  are  principally  related  to  obtaining  accurate  values  Jhe 

attenuation  of  the  amplitudes  of  Lg  waves,  so  that  mb(Lg)  can  be  estimated  to  closer 
tifn  onf-tLth  of  a  magnitude  unit.  Another  problem  is  to  obtain  data  or  e^losions 
of  announced  yield  for  diverse  areas  of  the  world,  to  test  the  universality  of  t 
yield  calibration  relations  for  all  continental  regions.  (continued) 
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19.  Abscract  (continued) 

The  calibration  equations  are: 

1)  for  water-saturated  rock 

mb(Lg)  =  3.943  +  1.124  log  Y  -  0.0829  (log  Y)^ 

with  a  standard  deviation  of  an  individual  value  from  the  curve  of  0.05,  or,  if  a 
linear  solution  is  preferred, 

mb(Lg)  =  4.307  ±  0.067  +  (0.765  +  0.027)  log  Y  for  5.2  <  mb  <  6.7 

2)  for  unsaturated  material 

mb(Lg)  =  3.869  +  1.110  log  Y  -  0.14  ii  (log  Y)^ 

with  a  standard  deviation  of  an  individual  value  from  the  curve  of  0.09,  or,  if  a 
linear  relation  is  preferred, 

mb(hg)  =  3.965  ±  0.049  +  (0.833  ±  0.848)  log  Y  for  4.0  £  mb  <  5.4 

When  these  equations  are  applied  to  seven  explosions  of  announced  yield  in  western 
United  States,  eastern  United  States,  French  Sahara  and  Shagan  River,  U.S.S.R.,  the 
yields  estimated  from  mb(Lg)  never  differ  from  the  announced  yields  by  more  than  36%. 

The  estimated  mb(P)  bias  between  NTS  and  eastern  North  America,  based  on  mb(Lg) 
data,  is  0.31  t  0.02  magnitude  units. 

Further  research  will  concentrate  on  estimating  yield  values  of  Soviet  explosions, 
and  of  estimating  the  mb(P)  bias  between  the  site  of  those  explosions  and  NTS. 

A  new  technique  which  combines  frequency -domain  Wiener  filtering  and  modal  isolation 
is  developed  for  calculating  interstation  phase  velocities,  group  velocities,  and 
attenuation  coefficients  of  seismic  surface  waves.  Frequency-domain  Wiener  filtering 
is  more  effective  than  time-domain  Wiener  filtering  for  these  determinations  because 
it  uses  a  small  window-lag  and  produces  a  smoother  interstation  Green's  function.  This 
leads  to  greater  accuracy  and  stability  when  noise-contaminated  data  are  analyzed. 

To  more  effectively  eliminate  the  effects  of  modal  interference,  phase-matched 
filtering  or  time-variable  filtering  can  be  applied  before  Wiener  filtering  to  isolate 
one  particular  mode  for  each  of  two  stations.  Synthetic  seismograms  contaminated 
by  noise  as  well  as  real  data  are  used  to  compare  this  technique  to  other  methods;  it 
is  found  to  be  superior  in  all  cases  where  the  wave  forms  are  contaminated  by  noise. 

The  application  of  Wiener  filtering  and  modal  isolation  to  surface  wave  data  in 
stable  continental  regions  produces  reliable  attenuation  coefficient  values  to  periods 
longer  than  those  which  have  previously  been  determined.  Inversion  of  those  data  leads 
to  lower  values  of  Qg  in  the  lower  crust  and  upper  mantle  than  those  obtain  in 
earlier  studies. 
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TECHNICAL  SUMMARY 


The  objectives  of  the  Lg  studies  of  underground  explosions  are  to 
obtain  calibration  curves  or  equations  of  m^(Lg)  versus  the  logarithm  of 
the  explosion  yield  and  to  obtain  estimates  of  the  mjj(p)  bias  between 
test  sites,  such  as  NTS  and  Shagan  Hiver,  by  means  of  ni^(Lg)  values. 
The  first  objective  will  provide  a  means  of  estimating  directly  the 
explosion  yield  from  measured  Lg-wave  amplitudes.  The  second  objective 
will  provide  the  ni^(P)  bias  information  that  is  needed  to  estimate 
yields  from  m^(P)  values  at  places  other  than  NTS. 

Technical  problems  are  principally  related  to  obtaining  accurate 
values  for  the  attenuation  of  the  amplitudes  of  Lg  waves,  so  that  m^(Lg) 
can  be  estimated  to  closer  than  one- tenth  of  a  magnitude  unit.  Another 
problem  is  to  obtain  data  for  explosions  of  Eunnounced  yield  for  diverse 
areas  of  the  world,  to  test  the  universality  of  the  yield  calibration 
relations  for  all  continental  regions. 

The  calibration  equations  are: 

1)  for  water-saturated  rock 

iBt,(Lg)  =  3.943  +  1.124  log  Y  -  0.0829  (log  Y)^ 
with  a  standard  deviation  of  an  individual  value  frcmi  the  curve  of  0.05, 
or,  if  a  linear  solution  is  preferred, 

»l,(Lg)  =  4.307  ±  0.067  +  (0.765  ±  0.027)  log  Y  for  5.2  i  m^  i  6.7 

2)  for  unsaturated  material 

mjj(Lg)  =  3.869  +  1.110  log  Y  -  0.146  (log  Y)^ 
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with  a  standard  deviation  of  an  individual  value  from  the  curve  of  0.09, 
or,  if  a  linear  relation  is  preferred, 

\(Lg)  =  3.965  ±  0.049  +  (0.833  ±  0.848)  log  Y  for  4.0  i  nij,  ^  5.4 

When  these  equations  are  applied  to  seven  explosions  of  announced 
yield  in  western  United  States,  eastern  United  States,  French  Sahara  and 
Shagan  River,  U.S.S.R.,  the  yields  estinated  from  m^(Lg)  never  differ 
from  the  announced  yield  by  more  than  36$. 

The  estimated  bias  between  NTS  and  eastern  North  America, 
based  on  m^(Lg)  data,  is  0.31  ±  0.02  magnitude  units. 

Further  research  will  concentrate  on  estimating  yield  values  of 
Soviet  explosions,  and  of  estimating  the  bias  between  the  site  of 
those  explosions  and  NTS. 

A  new  technique  which  combines  frequency-domain  Wiener  filtering 
and  modal  isolation  is  developed  for  calculating  Interstatlon  phase 
velocities,  group  velocities,  and  attenuation  coefficients  of  seismic 
surface  waves.  Frequency-domain  Wiener  filtering  is  more  effective  than 
time-domain  Wiener  filtering  for  these  determinations  because  it  uses  a 
small  window-lag  and  produces  a  smoother  interstation  Green's  function. 
This  leads  to  greater  accuracy  eind  stability  when  noise-contaminated 
data  are  analyzed.  To  more  effectively  eliminate  the  effects  of  modal 
interference,  phase-matched  filtering  or  time-variable  filtering  can  be 
applied  before  Wiener  filtering  to  Isolate  one  particular  mode  for  each 
of  two  stations.  Synthetic  seismograms  contaminated  by  noise  as  well  as 
real  data  are  used  to  compare  this  technique  to  other  methods;  it  is 
found  to  be  superior  in  all  cases  where  the  wave  forms  are  contaminated 
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by  noise 


The  application  of  Wiener  filtering  and  modal  isolation  to  surface 

wave  data  in  stable  continental  regions  produces  reliable  attenuation 

coefficient  values  to  periods  longer  than  those  which  have  previously 

been  determined.  Inversion  of  those  data  leads  to  lower  values  of  Q-  in 

P 

the  lower  crust  and  upper  mantle  than  those  obtained  in  earlier  studies. 
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YIELD  ESTIMATES  OF  NEVADA  TEST  SITE 
EXPLOSIONS  OBTAINED  FROM  SEISMIC  Lg  WAVES 
Otto  W.  Nuttll 


of  Earth  and 


Saint  Louis  University.  St.  Louis,  Missouri 


ABSTRACT 


A  methodology  is  presented  for  determining  the  yield  of  underground 
nuclear  explosions  from  Lg-wave  amplitudes.  The  methodology  is  applied 
to  Nevada  Test  Site  explosions,  for  which  the  data  from  short-period, 
vertical  component  analog  seismographs  at  three  stations  are  used  to 
develop  calibration  curves  for  unsaturated  material  and  water-saturated 
rock  source  conditions.  The  latter  curves  are  found  to  provide  reason¬ 
ably  accurate  estimates  of  the  yields  of  explosions  in  other  areas  of 
the  United  States  and  in  the  French  Sahara,  suggesting  that  they  may  be 
applicable  to  all  continental  areas.  If  so,  they  also  can  provide  an 
estimate  of  the  bias  of  nij^(P)  magnitudes  between  different  continental 
sites.  For  example,  the  Lg  data  from  NTS  explosions  indicate  a  0*31  ± 
0.02  magnitude  unit  bias  between  NTS  and  eastern  North  America,  similar 
to  the  approximately  0.33  unit  bias  found  between  western  and  eastern 
North  America  previously  by  use  of  earthquake  data. 


INTRODUCTION 


The  seismic  Lg  wave  is  one  of  a  number  of  regional  phases,  includ¬ 
ing  Pn,  Pg  and  Sn,  that  propagate  in  the  continental  lithosphere. 
Because  the  anelastic  attenuation  of  1-sec  period  Lg  waves  is  small  in 
shield  and  geologically  old  stable  regions,  Lg-wave  amplitudes  provide  a 
useful  tool  for  estimating  1-sec  period  magnitudes,  such  as  m^^,  for 
small  earthquakes  and  explosions  fNuttli.  1973].  Furthermore,  inasmuch 
as  Lg  consists  of  a  superposition  of  many  higher-mode  surface  waves  of 
group  velocity  near  3.5  km/sec,  its  radiation  is  more  isotropic  than 
that  of  P  and  S  waves.  This  feature  adds  to  its  usefulness  as  a  magni¬ 
tude  estimator  for  small  events,  because  full  azimuthal  coverage  is  not 
essential  and  thus  reliable  magnitude  determinations  can  be  made  from 
the  data  of  only  a  few  stations. 

For  routine  magnitude  estimates,  given  to  the  nearest  one- tenth  of 
a  unit,  reduction  of  amplitude  for  the  effect  of  epicentral  distance  can 
be  made  by  using  an  average  regional  value  of  the  coefficient  of  anelas¬ 
tic  attenuation.  As  examples,  for  eastern  North  America  the  coefficient 
of  anelastic  attenuation  of  1-sec  period  Lg  waves  may  be  taken  as  6  x 
10“^  km”^  rNuttll.  1973]  and  for  Iran  as  4.5  x  10”^  km"^  rNuttll.  198O]. 
However,  for  typical  versus  explosion  yield  relations,  e.g.  Bolt 
[1976],  Dahlman  and  Israelson  [1977],  Murphy.  [I98I],  a  one-tenth  unit 
error  in  corresponds  to  about  a  33%  error  in  explosion  yield.  There 
are  several  possible  approaches  for  reducing  the  portion  of  this  error 
that  results  from  Inexact  reduction  for  attenuation.  One  is  to  use  a 
regional  average  value  for  the  coefficient  of  anelastic  attenuation  and 
then  determine  source-to-station  corrections  for  each  station  with 


respect  to  a  given  source,  by  minimizing  the  sum  of  the  squares  of  the 
magnitude  residuals.  This  is  the  procedure  most  often  applied  for  mjj(p) 
and  Mg  determinations.  A  second  possible  approach  is  to  estimate  the 
value  of  the  coefficient  of  anelastic  attenuation  for  each  source- to- 
statlon  path,  as  a  function  of  wave  frequency.  In  this  paper  the  second 
approach  initially  is  used,  and  then  the  results  are  further  refined  by 
additional  application  of  the  first  approach. 

EXPERIMENTAL  DETERMINATION  OF 
ANELASTIC  ATTENUATION 

Data  5ei. 

The  Nevada  Test  Site  (NTS)  is  located  in  southern  Nevada  (approxi¬ 
mately  37°N,  116*^/).  Springer  and  Kinnaman  [19713  discuss  the  geography 
and  surficial  geology  of  the  test  site,  which  lies  within  the  young 
Basin  and  Range  province  of  the  western  United  States.  Values  of  the 
coefficient  of  anelastic  attenuation  of  1-sec  period  Lg  waves  in  the 
region  near  NTS  vary  from  4.7  to  6.0  x  10”^  km”^  F Singh  and  Herrmann. 
1983]f  based  on  data  obtained  from  the  World-Wide  Standardized  Seismo¬ 
graph  Network  (WWSSN)  and  the  Long  Range  Seismic  Measurements  (LRSM) 
network. 

Film  copies  of  short-period,  vertical-component  (SPZ)  WWSSN  seismo¬ 
grams  were  used  exclusively  in  this  study.  For  future  studies  these  data 
might  be  supplemented  by  those  of  other  networks,  such  as  the  LRSM  and 
the  University  of  California.  Lawrence  Livermore  National  Laboratory 
network. 

The  Lg-wave  amplitude  decrease  due  to  anelastic  attenuation  is 


proportional  to  exp(-)22k)>  where  ^  is  epicentral  distance  and  X  is  the 
coefficient  of  anelastic  attenuation,  related  to  the  quality  factor  Q  by 


y(f)  =  wf/U(f)  .  0(f)  (1) 

where  f  is  wave  frequency  and  U  is  the  Lg  group  velocity.  Assuming 
approximate  values  of  1  Hz  for  f,  3.5  km/ sec  for  U  (1  Hz)  and  150  for  Q 
(1  Hz)  =  0^.  it  follows  that  X  (1  Hz)  =  0.006  km""'.  At  ^  =  500  km.  a 
10?  error  in  would  correspond  to  a  21?  error  in  A(Lg),  where  A(Lg)  is 
the  Lg-wave  amplitude,  which  would  result  in  an  error  of  0.08  units  in 
the  1-Hz  Lg  magnitude.  At  ^  =  750  km,  a  10?  error  in 
translates  into  a  36?  error  in  A(Lg),  equivalent  to  a  0.13  unit  error  in 
in^(Lg)  or  about  a  50?  error  in  explosior  yield.  By  such  reasoning  data 
from  stations  beyond  about  750  km  from  NTS  were  not  used  in  this  study. 
This  restricted  the  selection  of  tfWSSN  stations  to  BKS  (Berkeley,  Cali¬ 
fornia),  DUG  (Dugway.  Utah),  GSC  (Goldstone,  California),  and  TUC  (Tuc¬ 
son,  Arizona).  Inspection  of  the  seismograms  for  GSC  revealed  that  they 
usually  were  off-scale  for  NTS  explosions.  This  also  occurred  at  TUC 
for  yields  of  100  kt  (kilotons)  and  greater  and  at  DUG  for  yields  of  50 
kt  and  greater,  except  when  the  magnification  of  the  seismograph  at  DOG 
was  reduced  prior  to  large  explosions.  Because  in  general  the  GSC 
seismograms  were  not  usable,  only  data  from  BKS,  DUG  and  TUC  were 
analyzed. 

Si  Deteriiiine4  Irom  Lg.  £.0<ia« 

The  frequency  dependence  of  Q  for  Lg  waves,  for  the  purposes  of 
this  study,  was  assumed  to  satisfy  the  relation  [MltQhsLLi  1980] 

Q(f)  =  Qj,  f^  (2) 
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where  ^  is  constant  over  the  range  of  frequencies  of  interest,  from 
about  2  to  0.1  Hz.  The  value  of  ^  and  the  first  approximation  to  the 


I 
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It 


value  of  for  paths  from  NTS  to  BKS,  DUG  and  TUC  are  obtained  by 
application  of  the  coda  Q  method  F Akl  and  Chouet.  1975]  as  adapted  by 
Herrmann  [1980].  By  this  method  the  coda  of  Q,  considered  to  be  pro¬ 
duced  by  scattered  waves,  ideally  will  exhibit  a  continuous  decrease  in 
wave  amplitude  and  wave  frequency  with  increasing  time  after  the  arrival 
of  the  Lg  wavetrain.  Measurements  of  the  wave  frequency  as  a  function 
of  time  after  the  origin  (lapse  time)  can  be  fitted  by  theoretical 
curves,  which  have  and  ^  as  parameters.  Applying  this  technique  to 
earthquake  data,  Singh  and  Herrmann  [1983]  obtained  values  of  150  to 
190  and  ^  values  of  0.5  to  0.6  in  the  area  near  NTS. 


According  to  theory,  the  frequency  of  the  coda  waves  should 

f 

decrease  as  the  lapse  time  increases  FAkl  aRi  Chouet.  1975;  Herrmann. 
1980].  In  actual  practice  there  are  two  types  of  observed  departure 
from  the  ideal  case.  The  first  type  is  caused  by  the  arrival  of 
^  fundamental-mode  Rayleigh  waves,  of  lower  frequency  eujd  often  larger 

amplitude  than  the  coda  waves,  with  group  velocities  between  about  2.6 
and  1.4  km/sec.  These  fundamental-mode  waves,  which  are  strongly  excited 

I 

by  shallow-depth  sources  such  as  explosions,  are  superposed  on  a  portion 
of  the  coda  in  the  time  window  in  which  they  arrive.  The  second  type  of 

departure  from  the  ideal  case  is  the  observed  tendency  of  the  coda  waves 

f 

to  be  of  constant  frequency  at  lapse  times  exceeding  a  value  t',  where 
t*  Increases  as  the  magnitude  of  the  earthquake  or  explosion  increases. 
The  dependence  of  t*  on  magnitude  suggests  that  the  phenomenon  in  some 

k 


way  is  related  to  the  source. 


Figures  1.  2  and  3  show  portions  of  the  BKS  seismograms,  SPZ  com¬ 
ponent,  for  the  NTS  explosions  STANLEY,  NASH  and  ALMENDRO,  respectively. 
Lapse  times  are  indicated  at  the  minute  marks.  Lower-case  letters  are 
used  to  indicate  wavelets  in  the  coda  that  are  used  for  estimating 
and  In  practice,  the  film  copy  of  the  seismogram  was  enlarged  to  6 
mm/sec  in  order  that  the  wave  periods  could  be  estimated  to  a  hundredth 
of  a  second.  Table  1  contains  the  lapse  times  and  wave  frequencies  of 
the  wavelets  to  be  used  in  estimating  the  values  of  and 

Fundamental-mode  wavelets  in  Figures  1  through  3  are  Indicated  by 
upper-case  letters.  Their  lapse  times  and  wave  frequencies  can  be  found 
in  Table  1 .  The  near-constant  frequency  wavelets  are  denoted  by 
numbers.  Their  lapse  times  and  wave  frequencies  also  are  included  in 
Table  1 .  Theory  indicates  that  the  fundamental-mode  data  should  not  be 
used  and  the  results  of  this  study  suggest,  as  will  be  shown  later,  that 
the  constant- frequency  data  should  not  be  used  in  the  curve-fitting  pro¬ 
cess  to  determine  Q  and  6  values. 

o  ' 

Figure  4  is  a  plot  of  the  lapse  times  and  wavelet  frequencies  as 
given  in  Table  1.  The  solid-line  cui’ve  is  the  author's  fit  of  a 
theoretical  curve  r Herrmann.  198O]  to  the  x  data  points  of  Table  1  and 
of  similar  data  for  15  additional  explosions.  The  parameters  of  this 
curve  initially  were  found  to  be  =  143.  ^  =  0.6  and  by  further 
refinement,  as  described  later,  were  found  to  be  s  139,  ^  :  0.6.  The 
dashed-line  curve  is  the  least-squares  fit  of  another,  but  similar,  set 
of  lapse  time  versus  wave  frequency  data  obtained  by  Rondout  r Peseokls 
and  Pomeroy.  1984;  Rondout  Associates.  1984],  with  =  225  and  ^ 

=  0.2.  This  curve  better  fits  the  constant-frequency  portion  of  the 


data  of  Figure  4,  as  indicated  by  o's,  than  the  lower  frequency  data 
points  indicated  by  x's.  The  fundastental-mode  data  points,  indicated  by 
•«-'s,  fall  well  below  the  curves. 

Singh  and  Herrmann  [1983]  gave  a  map  of  ^  values  for  the  United 
States,  as  determined  by  applying  the  coda  Q  r Herrmann.  198O]  method  to 
earthquake  data.  In  the  southern  portions  of  California  and  the  Basin 
and  Range  province  their  ^  values  lie  between  0.4  and  0.6.  In  the 
remainder  of  the  United  States  they  generally  are  between  0.1  and  0.3. 

Figure  5  shows  wave  frequency  versus  lapse  time  curves  for  BKS,  DUG 
and  TUC  for  explosive  sources  at  NTS.  A  fit  of  Herrmann* a  [198O]  curves 
to  the  data  gave  =  143  for  BKS,  150  for  DOC  and  I6O  for  TUC,  along 
with  ^  =  0.6  for  all  three  stations.  Further  minor  refinement,  as  dis¬ 
cussed  later,  resulted  in  the  values  corresponding  to  the  solid-line 
curves  in  Figure  5  (Q^  =  139  at  BKS,  155  at  DUG  and  162  at  TUC).  The 
dashed-line  curves  are  drawn  for  values  approximately  10)^  larger  and 
smaller  than  those  represented  by  the  solid-line  curves.  With  a  few 
exceptions,  the  dashed-line  curves  bound  the  data  points. 

The  assumed  exponential  dependence  of  Q  on  frequency,  as  given  by 
eq.  (2),  is  not  unique.  Rondout  Associates,  Inc.  [1984]  also  con¬ 
sidered  a  liner  relation  of  the  form 

Q(f)  =  +  c(  f 

They  found  for  BKS  that  their  data,  of  the  type  shown  in  Figure  4,  could 
be  fitted  equally  well  by  an  exponential  relation,  with  =  225  and  $  = 
0.2,  and  by  a  linear  relation,  with  &  175  and  <(  =  50. 


The  mjj(Lg)  scale,  as  originally  defined  FNuttll.  1973],  makes  use 

of  Lg-wave  amplitudes  at  periods  near  1  sec.  Therefore  an  accurate 

value  of  m^(Lg)  requires  an  accurate  value  of  Q^.  As  can  be  seen  from 

the  results  presented  above,  the  frequency  versus  lapse  time  data 

appear,  by  themselves,  to  be  inadequate  for  obtaining  a  unique  value  of 

Q  .  Therefore  additional  information  must  be  used  to  obtain  Q  values  by 
o 

independent  means.  Attempts  to  do  so  are  discussed  in  the  following 
three  sections. 

Measurements  jq£  Jissax.  s£.  L&  With  Pla..tano.e. 

The  first  of  these  methods  makes  use  of  the  amplitude  attenuation 
of  Lg  waves  as  a  function  of  epicentral  distance  for  waves  of  a  given 
frequency.  The  amplitudes  for  waves  of  a  selected  frequency  are  fitted 
by  a  theoretical  curve  which  takes  account  of  geometric  spreading, 
dispersion  and  anelastlc  attenuation,  and  has  y  as  a  parameter.  The 
average  value  of  in  the  area  containing  the  epicenter  and  stations  Is 
obtained  when  1-Hz  waves  are  used.  By  using  the  same  procedure  for 
waves  of  different  frequencies  the  dependence  of  Q  on  f.  or  the  value  of 
can  be  obtained.  The  method  works  well  when  is  approximately  con¬ 
stant  over  the  area,  as  It  Is  for  eastern  North  America  rNuttll.  1973], 
However,  for  the  NTS  region  the  amplitude  veCLues  were  found  to  scatter 
appreciably  due  to  the  variability  of  over  the  area,  so  that  one 

could  not  choose  between  the  two  estimates  obtained  from  the  coda 

o 

data. 

Meaattrcawita  SiL  \(£)  >  Bi,(Lg)  filaa  Xfill  Earthquakes. 


The  second  method  makes  use  of  earthquakes  that  occurred  near  NTS 


Dermenglan  et  al-  [1984]  provided  a  list  of  some  such  earthquakes.  If 
the  attenuation  of  P-wave  amplitudes  In  the  lithosphere  and  astheno~ 
sphere  below  NTS  Is  the  same  as  It  Is  beneath  eastern  North  America, 
earthquakes  or  explosions  of  the  same  strength  will  have  the  same  ni^(P) 
values.  However,  the  observed  mj^(p)  values  at  NTS  are  found  to  be 
smaller.  The  difference  In  values  Is  called  the  nijj(P)  bias  of  one 
region  with  respect  to  the  other.  Chung  and  Bernreuter  [1981],  using 
earthquakes,  found  the  bias  between  western  and  eastern  North  America  to 
be  approximately  0.33  body-wave  magnitude  units,  which  they  explained  as 
being  caused  by  Q  differences  In  the  asthenosphere. 

The  m^(Lg)  magnitude  scale  was  constructed  so  that  m^(P)  equals 
m^(Lg)  for  eastern  North  American  earthquakes  [Nuttll.  1973].  Lg-wave 
amplitudes  are  not  expected  to  be  affected  by  the  asthenosphere,  as  the 
wave  paths  are  confined  to  the  crust.  Also,  for  shallow  source  depths 
they  are  little  affected  by  variations  in  continental  crustal  structure 
r Camplllo  ^. ,  1984].  Therefore  m^(Lg)  Is  not  expected  to  show  bias, 
i.e.  for  the  same  strength  event  in  eastern  and  western  North  America 
the  m^(Lg)  values  should  be  the  same. 

Tedclng  into  account  the  finding  of  Chung  .aAd.  Bernreuter  [1981], 
m^(Lg)  ought  to  be  about  0.3  to  0.4  units  larger  than  mj^(P)  values  for 
earthquedces  near  NTS.  Table  2  shows  ni^(P)  values  for  seven  small  to 
moderate  Nevada  earthquakes.  For  most  of  them  the  m^(P)  value  is  based 
upon  the  amplitude  data  of  a  single  station,  making  them  less  reliable 
than  desirable.  Also  shown  In  Table  2  are  in^(Lg)  values  calculated 
assuming  both  the  Q^,  ^  values  of  Peseckls  ^Qd.  Pomerov  [1984]  and  those 
found  in  the  present  study.  The  former  set  of  values,  In  general. 


results  in  m^^CLg)  estimates  that  are  O.^l  to  0.5  units  smaller  than  the 
latter  set  of  estimates.  Furthermore,  the  Peseckis  and  Pomerov  [1984] 
values  of  Q^.  ^  in  all  cases  lead  to  m^(Lg)  estimates  that  are  smaller 
than  mjj(P),  rather  than  being  larger,  as  is  expected  because  of  the  P- 
wave  magnitude  bias.  For  four  of  the  seven  earthquakes  the  Q^,  i>  values 
of  the  present  study  result  in  the  mj^(Lg)  estimates  being  larger  than 
mj^CP),  as  expected.  For  two  of  the  earthquakes  the  opposite  is 
observed,  and  for  the  remaining  one  the  mjj(p)  and  mjj(Lg)  estimates  are 
essentially  the  same.  Although  not  as  conclusive  as  desirable,  the 
results  of  this  comparison  suggest  that  the  values  of  Peseckis  and 
Pomerov  [1984]  might  be  too  large  and  that  the  estimates  of  the 
present  study  are  closer  to  being  correct. 


Response  Spectrum 


Method. 


The  third  method  is  more  sensitive  to  ^  than  to  values. 
Although  la  more  important  than  $  for  the  purpose  of  this  paper, 
namely  to  determine  m^^CLg),  when  Q  l3  assumed  to  satisfy  an  exponential 
dependence  on  f  there  is  some  trad«j-off  between  and  ^  values.  The 
trade-off,  which  results  because  the  frequency  versus  lapse  time  data 
can  be  measured  over  only  a  limited  range  of  frequencies,  causes  the 
estimate  to  decrease  as  the  estimated  value  of  ^  Increases.  Therefore 
an  Independent  knowledge  of  ^  will  narrow  the  possible  range  of 
values.  For  example,  using  ^  s  0.6  in  Figure  5  restricts  the  values 
to  a  much  smaller  range  than  if  ^  were  also  permitted  to  vary. 


The  third  method  makes  use  of  response  spectra  that  are  calculated 
from  strong-motion  accelerograms  of  southern  California  earthquakes 
recorded  in  southern  California,  not  too  distant  fron  NTS.  Response 
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spectra  have  the  property  that  the  limiting  value  of  spectral  accelera¬ 
tion  at  high  frequencies  equals  the  maximum  observed  acceleration,  as 
read  directly  from  the  accelerogram  (see,  e.g.  Newmark  and  Hall 
[1982]).  Figure  6  shows  a  set  of  response  spectra,  for  5%  damping, 
extrapolated  from  an  eplcentral  distance  of  21 1  km  to  a  distance  of  50 
km  for  the  April  9,  1968  (02^29®  U.T. )  southern  California  earthquake 
(33.2°N,  116.1°W)  of  r  6.4.  Values  of  ^  between  0.4  and  0.9  were 

assumed,  along  with  a  of  150.  Lines  of  constant  spectral  accelera- 

—  S 

tlon,  from  10  g  to  10  g,  where  g  is  the  acceleration  of  gravity,  are 
plotted  in  the  figure.  The  best  value  of  ^  is  that  which  causes  the 
spectrum  at  short  periods  to  be  parallel  to  the  curves  of  constant 
acceleration,  and  to  have  a  limiting  value  of  spectral  acceleration 
numerically  equal  to  the  maximum  acceleration  as  read  directly  from  the 
accelerogram.  For  an  =  6.4  California  earthquake,  the  curves  of 
J-flYnsJC  ajai.  Sastca  C198I]  and  Campbell  [1981]  give  an  expected  a^^^  of 
0.04  g  to  0.05  g  at  50  km  distance.  In  Figure  6  the  ^  =  0.6  curve 
corresponds  to  an  a^^^  value  of  0.043  g.  The  calculated  spectra  for 
values  of  ^  less  than  0.6  curve  up  at  the  short  periods.  Indicating  that 

y  Is  too  large  and  ^  too  small.  For  $  =  0.7  the  spectral  curve  of  Fig¬ 

ure  6  yields  a^^  =  o.02  g,  which  fr<»i  Jovner  and  Boore  [1984]  and  Camo- 
Jiail  [1984]  corresponds  to  an  earthquake  of  s  5.O,  not  6.4.  For 

values  of  ^  larger  than  0.7  the  short-period  portion  of  the  spectra 
shown  in  the  figure  have  a  slope  greater  than  unity,  indicating  that  the 
assumed  y  value  Is  too  small  and  the  assumed  ^  value  too  large. 

Similar  analysis  was  used  for  eleven  other  response  spectra  of  the 
April  9,  1968  southern  California  earthquake  and  the  February  9,  197I 
San  Fernando,  California  earthquake.  In  all  eases  the  best  ^  values 
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were  found  to  lie  between  0.5  and  0.7,  with  a  median  of  0.6.  The  calcu¬ 
lations  were  repeated  for  assumed  values  of  200  and  250,  and  these 
particular  results  observed  to  be  relatively  insensitive  to  the  assumed 
Qq  value.  Therefore  the  appropriate  $  value  for  southern  California 
appears  to  be  0.6,  similar  to  that  obtained  in  the  present  study  for 
paths  from  NTS  to  BKS,  DOG  and  TOC.  Peseokls  and  Pomerov  [1984] 
obtained  ^  values  of  0.3  for  NTS  to  DOG  and  of  0.2  for  NTS  to  BKS  and 
TOC  when  they  fitted  their  coda-Q  data  by  an  exponential  model. 

mjj(Lg)  VERSOS  YIELD  CALIBRATION  CORVES 
FOR  NTS  EXPLOSIONS 

Springer  an!  Klnnaman  [1971.  1975]  published  yield  values  for  a 

number  of  O.S.  explosions,  inoludlng  some  of  those  not  at  the  NTS  site, 

along  with  hypocentral  coordinates  and  a  description  of  the  material  in 

which  the  explosion  was  fired.  They  noted  (pg.  1073):  *Yield8  are  the 

best  values  of  total  underground  energy  deposition  as  determined  by 

radiochemical  and  other  means.  They  are  expressed  in  kllotons  (kt), 

1 2 

where  one  klloton  is  defined  as  10  calories.  Absolute  accuracy  of  the 
radiochemical  measurements  is  not  known,  but  is  generally  assumed  to  be 
±.10  per  cent. 

For  the  purposes  of  this  study  the  NTS  explosion  media  were  divided 
into  two  principal  categories,  "water-saturated  rock"  and  "saturated 
material."  according  to  a  suggestion  of  R.G.  Cell  of  the  Lawrence  Liver¬ 
more  National  Laboratory  [personal  communication,  1985].  Springer  and 
ginnnman  [1971,  1975]  gave  the  depth  to  the  top  of  the  static  water 
table  and  the  depth  of  the  explosions,  and  in  some  instances  the  percent 
water  content  of  the  source  rock.  At  certain  points  of  NTS,  where  there 


is  a  parched  water  table,  there  may  be  a  significant  amount  of  water 
above  the  static  table. 

No  differentiation  between  explosion  sources  was  made  with  respect 
to  geographic  areas  of  NTS,  absence  or  presence  of  cratering,  or  rock 
type,  factors  that  are  known  to  affect  trij^CP)  estimates.  The  relatively 
small  number  of  explosions  of  announced  yield  at  NTS  precludes  such  dis¬ 
tinctions.  As  will  be  seen  later  In  this  section.  In  spite  of  Ignoring 
these  factors,  the  standard  deviation  of  a  single  mjj(Lg)  value  from  a 
mean  curve  of  m^^CLg)  versus  T  (logarithm  of  yield.  In  kllotons)  Is 
small,  about  0.05  to  0.06  body-wave  magnitude  units,  for  the  NTS  explo¬ 
sions  detonated  In  water-saturated  rock. 

The  third  largest  amplitude  In  the  time  window  corresponding  to 
group  velocities  of  3*6  to  3.2  km/sec  was  used  to  determine  m|j(Lg),  con¬ 
sistent  with  the  definition  of  the  m^j(Lg)  scale  (Nuttll.  1973].  The 
arrival  time  of  this  third  leu'gest  amplitude  eQso  was  measured,  for  the 
purpose  of  determining  the  average  group  velocity.  The  wave  period, 
when  It  could  be  measured,  always  was  found  to  lie  between  0.7  and  1.3 
sec,  consistent  with  the  definition  of  m^(Lg).  However,  wave  periods 
were  often  difficult  to  measure  because  of  large  amplitudes.  To  over¬ 
come  this  problem,  the  curves  of  Figure  5  were  used  to  obtain  the  wave 
frequency  corresponding  to  the  travel  time  (or  lapse  time)  of  Lg.  The 
values  of  Q  and  of  y  at  that  frequency  were  calculated  by  means  of  equa¬ 
tions  (2)  and  (1).  Next  the  Lg  amplitude,  eifter  correction  for  Instru¬ 
ment  magnification,  was  extrapolated  to  a  hypothetical  vzdue  at  a  dis¬ 
tance  of  10  km  by  means  of  the  formula,  adapted  from  that  on  pg.  358  of 
Ewing  et  al.  [1957] 
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(3) 


Ado  km)  =  A(^  [sln(^111.1)/alndO/111.1) 

expCyQ;;^  -  10)3 

where  AdO  km)  Is  the  hypothetical  Lg  amplitude  at  10  km  eplcentral  dis¬ 
tance  and  A(^  is  the  observed  amplitude  at  distance  ^  In  kilometers. 
The  term  (^10)^^^  is  the  correction  required  for  dispersion  when  ampli¬ 
tude  measurements  are  made  in  the  time  domair.  Camplllo  et  al.  [1984] 
showed  by  synthetic  seismograms  of  Lg  that  Uie  exponent  should  be  1/3, 
corresponding  to  an  Airy  phase,  rai,her  than  1/2,  corresponding  to  nor¬ 
mally  dispersed  surface  waves.  Nuttll  [1973]  came  to  a  similar  conclu¬ 
sion  from  observational  data  of  earthquakes.  The  term 
[sln(^111.1)/slnd0/111.1)]^^^  corrects  for  geometrical  spreading,  and 
the  term  exp[)'(^-  10)]  for  anelastic  attenuation. 

Nuttli  [1973]  set  the  reference  level  of  the  mjj(Lg)  scale  such  that 
an  m^(P)  s  5.0  earthquake  in  eastern  North  America  gave  a  hypothetical 
1-sec  period,  vertical-component  Lg  amplitude  of  110  micrometers  at  10 
km  eplcentral  distance.  From  this  it  follows  that 

nj,(Lg)  =  5.0  +  log^pCAdO  km)/110]  (4) 

where  AdO  km)  is  in  microueters  of  ground  motion. 

By  means  of  equations  (1)  through  (4),  m^(Lg)  values  were  calcu¬ 
lated  using  the  data  of  BKS,  DUG,  and  TOC  for  35  NTS  explosions  of 
announced  yield  and  for  SHOAL,  which  was  in  Nevada  but  off  the  test 
site.  Next  a  reestimation  procedure  was  used,  consisting  of  three 
steps. 

1 .  The  water- saturated  rock  and  unsaturated  material  data 
separately  were  fitted  by  least  squares  to  quadratic  curves.  A 


quadratic  curve,  rather  than  a  linear,  was  selected  to  account  for  the 
fact  that  the  slope  of  the  nijj(Lg)  versus  log  Y  curve  should  be  1.0  when 
the  corner  period  of  the  Lg  spectrum  is  much  less  than  1  sec  and  should 
be  0.5  when  the  corner  period  is  much  greater  than  1  sec. 


2.  Station  corrections  were  obtained  for  BKS,  DOG  and  TOC.  For  an 
individual  station  the  deviations  of  its  mjj(Lg)  values  with  respect  to 
the  avereige  of  the  mj^(Lg)  values  for  the  three  stations  were  calculated. 
Then  the  original  value  obtained  by  the  coda-Q  method  was  adjusted  so 
as  to  make  the  arithmetic  average  of  the  deviations  equal  to  zero.  The 
adjustments  were  small,  namely  from  1l»3  to  139  for  BKS,  150  to  155  for 
DUG  and  160  to  162  for  TUC.  This  indicates,  at  least,  a  high  precision,' 
or  relative  accuracy,  of  values  obtained  by  the  coda-Q  method.  The 
corrections  affect  the  individual  m^^CLg)  estimates  by  no  more  than  ±.0.04 
magnitude  units. 


3.  Finally,  m^(Lg)  values  were  recomputed  for  all  36  explosions. 
The  results  are  given  in  Table  3,  along  with  the  m^j(P)  values  published 
by  the  International  Seismological  Centre  (ISC)  and  the  yield  values 
published  by  Springer  and  Klnnaman  [1971.  1975]. 

The  results,  presented  in  Table  3,  along  with  the  quadratic  least- 
squares  curves,  are  plotted  in  Figures  7  and  8.  The  empirical  equation 
of  the  water-saturated-rock  calibration  curve  is 


i\(Lg)  =  3.943  +  1.124  log  Y  -  0.0829  (log  Y)^  (5) 

where  Y  is  in  kllotons  (kt).  The  slope  of  the  curve  at  Y  =  10  kt  is 
0.96,  at  Y  =  100  kt  is  0.79 »  and  at  Y  =  1000  kt  is  0.63.  As  noted  ear¬ 
lier.  the  slope  of  the  curve  is  expected  to  be  1 .0  for  small  yield 


explosions,  whose  spectral  corner  period  is  less  than  1  sec,  and  is 
expected  to  be  0.5  at  large  yields,  assuming  the  yield  is  linearly 
related  to  the  seismic  moment. 

The  anpirlcal  equation  of  the  unsaturated  material  calibration 
curve  is 


Di^(Lg)  =  3.869  +  1.110  log  Y  -  0.146  (log  Y)^  (6) 

where  the  explosions  of  10  kt  and  greater  were  given  twice  the  weight  of 
the  smaller  ones,  as  the  latter  had  small  Lg  amplitudes.  The  slope  of 
this  curve  is  1.11  at  Y  =  1  kt,  0.82  at  Y  r  lo  kt,  and  0.53  at  Y  =  100 
kt. 


The  standard  deviation  of  the  20  water-saturated  rock  data  points 
(PILE  DRIVER  and  SHOAL,  in  granite,  were  excluded)  from  the  curve  speci¬ 
fied  by  equation  (5)  is  0.05  magnitude  units.  The  largest  positive 
deviations  are  for  PILE  DRIVER  (+0.19),  SHOAL  (+0.12)  and  MINIATA 
(+0.15).  The  largest  negative  deviations  are  for  COMMODORE  (-0.12)  and 
MISSISSIPPI  (-0.07).  The  standard  deviation  of  the  13  unsaturated 
material  data  points  from  equation  (6)  is  0.09  magnitude  units.  The 
largest  positive  deviations  are  for  PALANQUIN  (+0.15)  and  SCHOONER 
(+0.12).  The  largest  negative  deviations  are  for  POMMARD  (-0.13)  and 
MERLIN  (-0.12). 

Empirical  formulas  for  determining  explosion  yield  by  use  of  either 
m^(P)  or  values  usually  are  given  in  a  linear  form,  e.g.  Dahlman  and 
Israelson  [1977,  pg.  269],  rather  than  the  quadratic  relation  assumed 
for  eqs.  (5)  and  (6).  If  a  llenar  relation  is  assumed,  the  values  from 
Table  3  give  the  alternate  equations: 
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water- saturated  nofiJi 


\(Lg)  =  1*.307  ±0.067  ♦  (0.765  ±  0.027)log  Y 

for  5.2  i  i  6.7  (5a) 

unsaturated  maJ^fiHlal 

“b(Lg)  =  3.965  ±  0.0M9  +  (0.833  ±  0.0J<8)log  Y 

for  it.O  i  i  5.4  (6a) 

where  the  quantities  following  the  ±  signs  are  standard  deviations  of 
the  Intercept  and  slope,  respectively.  Linear  fits  to  the  data  of 
water-saturated  rock  and  unsaturated  material  are  shown  In  Figures  9  and 
10,  respectively. 

The  deviation  of  mjj(Lg)  values  of  explosions  from  mean  curves,  such 
as  that  of  Figure  7,  can  be  used  to  study  the  effect  of  release  of  tec¬ 
tonic  strain  energy  that  accompanies  large  explosions.  Svkes  and 
Clufentes  [1984]  discussed  the  effect  of  this  phenomenon  on  vedues. 

Wallace  [1985]  noted  that  tectonic  strain  release  for  NTS  explosions 

appears  to  be  frequency  dependent.  His  review  of  the  literature  on  the 
subject  Indicated  that  there  Is  little  discernible  effect  on  the  ampli¬ 
tude  of  short  period  waves.  He  noted  that  there  was  a  strong  tectonic 
release  signature  for  longer  period  waves  (^  5  sec)  for  certain  NTS 
explosions.  He  determined  seismic  moments  from  SH  waves  for  21  NTS 
events.  Six  of  these  explosions  were  of  announced  yields.  Table  4 
lists  their  moments  and  deviations  of  m^(Lg)  from  the  mean  curve  of  Fig¬ 
ure  7.  There  Is  no  obvious  correlation  between  the  and  ^^(Lg) 
values  In  Table  4,  which  suggests  that  the  effect  of  tectonic  strain 
release  on  m^(Lg),  if  any.  Is  smaller  than  other  factors  responsible  for 
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the  scatter  in  mjj(Lg)  values. 


Use  of  values  of  and  ^  other  than  those  proposed  in  this  study 
will  lead  to  different  numerical  coefficients  in  equations  (5),  (5a), 
(6)  and  (6a).  As  long  as  the  sets  of  0^,  ^  values  for  the  various 
source-to-statlon  paths  are  consistent  among  themselves,  such  as  are 
those  proposed  by  Rondout  Associates.  Inc.  [198i»],  the  resulting  empir¬ 
ical  in^(Lg)  versus  yield  equations  will  give  reliable  estimates  of  yield 
for  NTS  explosions,  as  demonstrated  by  Rondout  Associates.  Inc.  [1984], 
Likewise  the  standard  deviation  of  the  individual  nij^(Lg)  values  from  the 
average  curve  will  be  small.  Absolute  accuracy  in  injj(Lg)  values  is 
required  only  when  a  calibration  curve  or  equation,  developed  from  data 
of  one  area,  is  applied  to  data  from  another  8u?ea. 

mjj(P)  BIAS  BETWEEN  NTS  AND  EASTERN 
NORTH  AMERICA 

If  the  Qq,  $  values  for  NTS  events  shown  in  Figure  5  are  correct, 
then  m^(Lg)  values  for  NTS  explosions  and  earthquakes  might  be  assumed 
to  be  the  same  as  for  similar  size  explosions  and  earthquakes  in  eastern 
North  America  (ENA),  l.e.  “b(Lg)jjYs  =  “b^^®^ENA*  However,  the  mjj(P) 
values,  obtained  from  teleseismlc  P-wave  amplitudes,  can  be  expected  to 
be  different.  Therefore,  the  ni|j(P)  bias  between  NTS  and  ENA  can  be 
expressed  as 

“b^^^ENA  "  “b^^^NTS  "  “b^^®^ENA  "  “b^^^NTS  “ 

“b^Lg)uTS  ■  “b^**^NTS 

because  “jj(Lg)gjjj^  was  defined  to  be  equal  to  Ob(P)ENA  and 

assumed  to  be  Independent  of  geographic  region.  Table  5,  together  with 


•  ^  >”•  fc'*  >.**•  k 
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Table  3,  contains  86  explosions  for  which  both  an  mjj(Lg)  and  an 
ISC)  are  available.  The  average  values  of  the  bias  of  NTS  with 

respect  to  ENA  is  found  to  be  O.31  ±  0.02  magnitude  units,  where  the 
0.02  is  the  standard  error  of  the  mean.  This  is  in  good  agreement  with 
the  value  of  approximately  0.33  found  previously  by  Chung  and  Bernreuter 
[1981]  using  earthquake  data. 

OTHER  D.S.,  FRENCH  SAHARA  AND  EAST  KAZAKH 
EXPLOSIONS  OF  ANNOUNCED  YIELD 

There  are  only  a  few  published  yield  values,  in  addition  to  those 
for  NTS,  for  underground  explosions  in  continental  regions.  They  £u*e 
SALMON  in  Mississippi.  GASBUGGY  in  New  Mexico  and  RULISON  in  Colorado 
r Springer  and  Kinnaman,  1971]t  RIO  BLANCA  in  Colorado  r Springer  and  Kin- 
naniAn,  1975],  RUBIS  and  SAPHIR  in  the  French  Sahara  r Marshall  aJl 
19791 >  and  four  at  widely  separated  points  in  the  U.S.S.R.  r Marshall  at 
al. .  1979].  If,  by  analogy  with  earthquakes,  m|^(Lg)  is  independent  of 
crustal  structure,  (l.e.  m^(Lg)  for  eastern  North  American  earthquakes 
equals  m^^CLg)  for  western  North  American  earthquakes),  the  mj^(Lg}  versus 
explosion  yield  relations  derived  from  NTS  data  might  be  expected  to 
apply  to  all  continental  regions.  To  test  this  hypothesis,  m^(Lg) 
values  are  calculated  for  the  events  listed  above,  except  for  three  of 
the  explosions  in  the  U.S.S.R.  for  which  no  seismograms  were  available 
to  the  author. 

Table  6  contains  the  results  of  applying  the  coda-Q  method 
r Herrmann.  198O]  to  the  non-NTS.  United  States  explosions.  In  spite  of 
the  relatively  large  range  in  values  and  eplcentral  distances  for  the 
SALMON  to  station  paths,  the  individual  m^(Lg)  values  do  not  differ  much 


from  the  average  value.  For  the  three  western  explosions  the  seismo¬ 
grams  of  nearby  stations,  GOL  and  ALQ,  were  off-scale.  For  each  of 
these  explosions  only  one  usable  WWSSN  station  could  be  found. 

No  information  concerning  water  content  of  the  salt  in  which  SALMON 
was  fired  is  given  by  jSprlhgfiC.  and  Klnnaman  [1971].  Accordingly,  in 
Table  6  yields  are  estimated  from  its  m^^CLg)  value  using  both  the  equa¬ 
tions  for  water-saturated  rock  and  for  unsaturated  material.  The  former 
gives  a  value  slightly  less  than  the  announced  yield  (11^),  whereas  the 
latter  gives  values  somewhat  greater,  or  28$.  The  other  three 

explosions  of  Table  6  can  be  considered  to  have  been  emplaced  in  water- 
saturated  rock.  For  all  these  events,  equation  (5a)  gives  smaller  yield 
estimates  than  (5).  The  greatest  difference  between  announced  and 
estimated  yield  is  for  RIO  BLANCA,  when  using  equation  (5),  for  which 
the  difference  is  2^  of  the  announced  value. 

The  similarity  between  the  announced  yields  and  the  m^(Lg)  yields, 
as  estimated  using  equations  (5)  or  (6),  suggests  that  the  yield  versus 
calibration  curves  for  NTS  can  be  applied  to  both  the  Rocky  Moun¬ 
tain  and  the  Gulf  Coast  areas. 

Assuming  m^(Lg)  to  be  equal  numerically  to  d^(P)  in  eastern  North 
America,  from  Table  6  the  magnitude  bias  between  the  sites  of  GASBUGGY, 
RULISON  and  RIO  BLANCA  with  respect  to  eastern  North  America  is  approxi¬ 
mately  0.6  magnitude  units,  similar  to  the  value  obtained  by  Raoof  and 
Nuttll  [1984]  for  portions  of  the  west  coast  of  South  America.  On  the 
other  hand,  the  magnitude  bias  between  the  site  of  SALMON  and  eastern 
North  America  can  be  seen  to  be  only  -0.06  magnitude  units,  or  near 
zero,  which  implies  that  the  excitation  of  1-Hz  Lg  waves  by  explosions 


'■j'  "r.'  ^ 


and  earthquakes  of  equal  o^^CP)  values  In  eastern  North  America  is  the 
same. 

Table  7  presents  estimates  of  the  yields  of  the  French  Sahara 
explosions  RUBIS  and  SAPHIR  and  of  one  cratered  explosion  at  East 
Kazakh,  U.S.S.R.  Equations  (b)  and  (5a)  also  give  reasonably  close 
estimates  of  the  announced  yields,  as  noted  earlier  for  non-NTS  explo¬ 
sions  in  the  United  States.  Of  the  six  estimated  yield  values  given  in 
Table  7,  the  differences  from  the  announced  yields  lie  between  1)%  and 
35%. 


The  magnitude  bias  of  the  Sahara  site  with  respect  to  eastern  North 
America  is  about  -0.1  to  -0.2  magnitude  units,  based  on  the  values  for 
RUBIS  and  SAPHIR.  The  negative  value  may  result  from  the  fact  that  the 
explosions  were  fired  in  granite.  Because  the  mj^(P)  values  of  the 
U.S.C.G.S,  and  I.S.C,  are  so  greatly  different  for  the  East  Kazakh 
explosion  (6.3  and  5.8,  respectively),  no  estimate  of  the  bias  between 
that  site  and  eastern  North  America  can  be  made  from  the  explosion  of 
January  15.  1965. 


DISCUSSION  AND  CONCLUSIONS 

A  methodology  for  estimating  yield  from  Lg-wave  amplitudes  has  been 
presented  for  the  case  when  the  data  from  only  a  few  stations  (three  or 
less)  are  available  for  a  region  of  rapidly  varying  Q.  These  are  the 
least  desirable  conditions  for  estimating  explosion  yields  frcmi  Lg 
amplitudes.  Better  estimates  of  in^(Lg)  and  of  yield  could  have  been 
obtained  If  the  data  of  more  stations  had  been  used,  or  if  the  region 
surrounding  the  test  site  were  one  of  higher  and  more  uniform  Q  values, 


such  as  eastern  North  America  or  central  Asia.  In  spite  of  this  the 
scatter  of  data  points  about  the  average  water-saturated  rock  curve  in 
Figure  7  is  small  (standard  deviation  equal  to  0.05  magnitude  units)  and 
may  be  more  the  result  of  variations  in  rock  type  in  the  source  region 
than  of  ramdom  scatter.  The  fact  that  the  methodology  results  in  rea¬ 
sonably  accurate  estimates  of  both  the  yield  and  the  P-wave  magnitude 
bias,  with  respect  to  eastern  North  America,  of  explosions  in  such 
diverse  areas  as  Colorado,  New  Mexico,  Mississippi  and  the  Sahara,  gives 
confidence  in  the  applicability  of  the  methodology  and  of  equations  (5) 
or  (5a)  to  other  continental  areas.  However,  the  number  of  available 
explosions  to  test  this  conclusion  was  small,  and  as  a  result  it  must  be 
considered  as  tentative. 

The  values  of  the  coefficients  of  log  Y  and  (log  Y)^  in  equation 
(5)  could  be  refined  by  use  of  additional  NTS  data,  namely  those  of 
explosions  of  unannounced  yield.  This  particularly  would  be  desirable 
for  explosions  of  1  kt  yield  and  less  in  water-saturated  rock,  because 
the  slope  of  equation  (5)  at  those  yields  exceeds  unity,  suggesting  that 
equation  (5)  may  overestimate  the  yields  of  small  explosions,  those  of 
yield  less  than  10  kt.  If  a  linear  relation  is  preferred,  such  as  eq. 
(5a),  the  slope  and  intercept  should  be  determined  separately  for  explo¬ 
sions  of  yield  less  than  and  greater  than  10  kt. 

If.  as  suggested  by  the  data  presented  in  this  paper,  equation  (5) 
or  (5a)  is  applicable  to  other  water- saturated  rock  continental  sites, 
then  the  use  of  additional  NTS  data  to  refine  the  calibration  relation 
can  be  Justified.  However,  if  the  coefficients  of  the  terms  in  equation 
(5)  vary  slightly  from  one  geographic  area  to  another,  then  the  results 
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obtained  in  this  paper  are  indicative  of  the  type  of  accuracy  that  might 
be  achieved  if  there  were  a  threshold  treaty  that  provided  for  only  a 
limited  number  of  calibration  explosions  and  a  limited  number  of  seismo¬ 
graph  stations. 

One  of  the  surprising  results  of  this  study  is  that  seismographs 
with  limited  dynamic  range,  such  as  the  WWSSN  short-period  instruments, 
can  provide  yield  estimates  from  explosions  as  small  as  1  kt  and  as 
large  as  1000  kt.  However,  at  the  1  kt  level  the  oberved  signal-to- 
nolse  ratio  for  Lg  was  small,  and  at  the  1000  kt  level  the  Lg  waves  fre¬ 
quently  were  off-scale  for  the  analog  seismograms.  The  use  of  digital 
seismographs  with  a  broader  dynamic  range  is  desirable  for  future  work. 

The  greatest  uncertainty  in  Lg  estimates  of  yield  results  from 
uncertainty  in  the  values  of  and  $  that  are  required  for  correcting 
the  observed  Lg  amplitudes  for  anelastic  attenuation.  In  this  paper  the 
coda-Q  method  as  presented  by  Herriaaha  [1980]  was  the  primary  means  for 
estimating  the  values  of  the  parameters.  Even  relatively  large  errors 
in  the  values  of  these  parameters  can  be  tolerated  if  a  different  cali¬ 
bration  equation  is  used  for  each  test  site  [Rondout  Associates.  Jqs..  , 
1984].  However,  if  a  universal  calibration  equation  is  desired,  and  if 
the  m^(Lg)  values  are  also  to  be  used  for  estimating  m^(P)  bias  between 
different  test  sites,  then  accurate  absolute  values  of  and  ^  are 
needed.  Independently  of  the  coda-Q  method,  they  can  be  obtained  by 
measuring  Lg  amplitudes  at  various  frequencies,  as  recorded  by  broad¬ 
band  Instruments,  along  a  profile  of  temporary  stations  located  between 
the  test  site  and  the  permanent  stations.  If,  for  some  reason,  such 
experiments  cannot  be  carried  out,  then  one  can  reaort  to  the  indepen- 


dent  methods  used  in  this  paper,  namely  the  use  of  strong-motion 
response  spectra  and  the  comparison  of  ni^(P)  and  m^(Lg)  for  earthquakes 
in  the  area  of  the  test  site. 

The  explosion  data  used  in  this  paper  give  a  value  for  the  ni^(p) 
bias  between  NTS  and  eastern  North  America  of  O.3I  "  0.02  magnitude 
units.  This  value  is  similar  to  that  obtained  by  Chung  and  Bernreuter 
fl98l]  from  earthquake  data,  who  obtadned  a  bias  of  approximately  0.33 
magnitude  units  between  western  and  eastern  North  America. 
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FIGURE  CAPTIONS 


Figure  1.  Copy  of  a  portion  of  the  3KS  short-period,  vertical-component 
seismogram  for  the  NTS  explosion  STANLEY  in  dry  tuff.  The  seismo¬ 
gram  copy  has  been  cut  into  strips  for  display  purposes,  so  as  to 
include  the  time  interval  of  1  to  6  min  after  the  explosion.  In 
the  portion  between  2  and  3  1/2  min  the  seismogram  trace  has  been 
enhanced  by  pencil  to  make  it  more  legible.  Lower-case  letters 
refer  to  wavelets  which  are  used  to  obtain  and  ^  from  the  plot 
of  wave  frequency  versus  lapse  time  (Figure  6).  Wavelets  marked  by 
upper-case  letters  are  believed  to  be  fundamental-mode  Rayleigh 
waves,  and  thus  are  not  used  to  obtain  and 

Figure  2.  Copy  of  a  portion  of  the  BKS  short-period,  vertical-component 
seismogram  for  the  NTS  explosion  NASH  in  dolomite.  The  figure 
includes  the  time  interval  of  1  to  7  min  after  the  explosion.  The 
numbers  1  through  3  at  the  lower  right  hand  side  refer  to  wavelets 
at  the  end  of  the  coda  whose  frequency  remains  nearly  constant  as 
the  lapse  time  increases.  They  are  not  used  to  obtain  and  ^ 
values.  See  the  caption  of  Figure  1  for  further  explanation  of 
symbols. 

Figure  3.  Copy  of  a  portion  of  the  BKS  short-period,  vertical-ccxnponent 
seismogram  for  NTS  explosion  ALMENDRO.  The  figure  includes  the 
time  interval  of  6  to  16  min  after  the  explosion.  The  wave  amplll- 
tudes  in  the  preceding  minutes  are  too  large  to  be  reproduced  con¬ 
veniently.  See  the  captions  of  Figures  1  and  2  for  further  expla¬ 
nation  of  symbols. 
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Figure  4.  Plot  of  the  wave  frequency  versus  lapse  time  data  for  the  Lg 
coda  of  the  BKS  seismograms  of  Figures  1  through  3.  The  x's 
correspond  to  ■'’e  lower-case  letters  in  these  figures,  the  +*s  to 
fundamental-mode  Rayleigh  waves  and  the  o's  to  constant  frequency 
wavelets  at  the  end  of  the  coda.  The  solid-line  curve  is  taken  to 
be  the  approximate  fit  to  the  data  (the  x's).  The  dashed-llne 
curve  is  that  obtained  by  Rondout  Associates.  Inc.  [1984]  for  NTS 
events  recorded  at  BKS. 

Figure  5.  Wave  frequency  versus  lapse  time  data  for  selected  NTS  explo¬ 
sions  recorded  at  DUG,  BKS  and  TUC.  The  data  originally  were  fit¬ 
ted  by  curves  of  ^  =  0.6  and  =  150,  143.  and  160,  respectively. 
The  solid-line  curves  plotted  in  the  figure  correspond  to  values 
after  station  corrections  were  made,  as  discussed  in  the  text.  The 
dashed-line  curves,  which  are  for  values  approximately  lOJ 
greater  than  and  lesr  than  the  mean  values,  bound  most  of  the  data 
points. 

Figure  6.  5%  damped  pseudovelooity  response  spectra  obtained  by  extra¬ 

polating  an  observed  spectrum  at  a  distance  of  211  km  to  50  km, 
assuming  =  150  and  various  values  of  ^  between  0.4  and  0.9.  The 
dashed  lines  are  lines  of  constant  response  spectrum  acceleration. 
The  e2u'thquake  occurred  in  southern  California  on  April  9,  1968, 
and  had  =  6.4  and  epioentral  coordinates  of  33.2°N,  116.1®W. 
The  coordinates  of  the  accelerograph  were  34°08'N,  118°07'W.  An 
upturn  of  the  spectra  at  short  periods,  such  as  for  ^  =  0.4  and 
0.5,  indicates  that  the  assumed  value  of  ^  is  too  small.  A  slope 
greater  than  unity  at  the  short  periods  Indicates  that  the  assumed 


value  of  ^  ia  too  large. 


Figure  7.  Jn^(Lg)  versus  logarithm  of  axplosion  yield  for  NTS  explosions 
of  announced  yield  in  water- saturated  rock.  The  curve  is  a 

second-degree  polynomial  obtained  by  least-squares  fit  to  the  data, 
and  is  given  by  equation  (5). 

Figure  8.  ai^j(Lg)  versus  logarithm  of  explosion  yield  for  NTS  explosions 
of  announced  yield  in  unsaturated  material.  The  curve  is  a 

second-degree  polynomial  obtained  by  weighted  least-squares  fit  to 
the  data,  and  is  given  by  equation  (6). 

Figure  9*  o^^CLg)  versus  logarithm  of  explosion  yield  for  NTS  explosions 
of  announced  yield  i  i  water-saturated  rock.  The  curve  is  a  linear 
least-squares  fit  to  the  data,  and  is  given  by  equation  (5a). 

Figure  10.  “jj(Lg)  versus  logarithm  of  explosion  yield  for  NTS  explo¬ 
sions  of  announced  yield  in  unsatuiated  materials.  The  curve  is  a 
linear  least-squares  fit  to  the  data  and  is  given  by  equation  (6a). 
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TABLE  2.  Magnitudes  of  Nevada  earthquakes  near  NTS 


Date 

Station 

mb(lg) 

(Rondout  t  values) 

mb(Lg) 

(Q„*y  values  of 
°this  study) 

*  * 

07/18/63 

3.9®" 

DUG 

3.33 

3.68 

07/20/63 

4.1^ 

BKS 

DUG 

3*69 

4.15 

4.31 

4.51 

11/16/63 

4.1^ 

BKS 

DUG 

3.08 

3-^9 

3.54 

3.90 

^  jii . 

08f2ll6k 

3.8^ 

DUG 

TUG 

3.^5 

3.03 

3.81 

3.74 

11/17/65 

3.?^ 

DUG 

TUG 

3.67 

3.02 

3.94 

3.82 

04/06/66 

4.1^ 

TUG 

3.48 

4.23 

^  — *  ■ 

fc  . 

08/16/66 

5.6^ 

BKS 

TUG 

5.03 

5.64 

5.77 

6.34 

^  Value  from  U.S.C.G.S. 
^  Value  from  I.S.C. 


TABLE  3*  Magnitudes  and  yields  of  NTS  explosions  of  announced  yield 


Date 

Name 

Rock  Type 

Stations  Used 

m^(P-ISC) 

Announced 

Yield^  (kt) 

06/27/62 

HAYMAKER 

dry  alluvium 

BKS 

5.42 

67 

1 0/05/62 

MISSISSIPPI 

tuff 

BKS 

5.82 

110 

09/13/63 

BILBY 

wet  tuff 

BKS 

6.09 

— 

235 

10/26/63 

SHOAL 

granite 

BKS, 

TUC 

5.19 

— 

12 

10/09/64 

PAR 

dry  alluvium 

BKS, 

DUG, 

TUC 

5.22 

4.8 

38 

11/05/64 

HANDCAR 

dolomite 

BKS, 

DUG, 

TUC 

5.03 

4.8 

12 

12/16/64 

PARROT 

dry  alluvium 

DUG, 

TUC 

3.91 

— 

1.2 

12/16/64 

MUDPACK 

dry  tuff 

BKS, 

TUC 

4.32 

.... 

2.7 

02/16/65 

MERLIN 

dry  alluvium 

BKS, 

DUG, 

TUG 

4.71 

10 

04/14/65 

PALANQUIN 

dry  rhyolite 

BKS, 

DUG, 

TUC 

4.66 

— 

4.3 

06/11/65 

PETREL 

dry  alluvium 

DUG, 

TUC 

4.02 

- - 

1.2 

02/24/66 

REX 

wet  tuff 

BKS, 

TUC 

5.28 

5.0 

16 

04/14/66 

DURYEA 

rhyolite 

BKS, 

DUG, 

TUC 

5.66 

5.4 

65 

05/05/66 

CYCLAMEN 

dry  alluvium 

BKS, 

DUG, 

TUC 

4.85 

4.4 

13 

05/ 06/ 66 

CHARTREUSE 

rhyolite 

BKS, 

DUG 

5.76 

5.4 

70 

05/27/66 

DISCUS  THROWER 

tuff 

BKS, 

TUG 

5.25 

5.0 

21 

06/02/66 

PILE  DRIVER 

granite 

BKS, 

DUG, 

TUC 

5.84 

5.6 

56 

06/30/66 

HALF  BEAK 

rhyolite 

BKS, 

TUC 

6.26 

6.1 

300 

12/20/66 

GREELEY 

wet  tuff 

BKS, 

DUG 

6.49 

6.3 

825 

05/20/67 

COMMODORE 

wet  tuff 

BKS 

6.04 

5.8 

250 

05/23/67 

SCOTCH 

wet  tuff 

BKS, 

DUG, 

TUC 

6.08 

5.7 

150 

05/26/67 

KNICKERBOCKER 

rhyolite 

BKS, 

DUG 

5.75 

5.4 

71 

09/21/67 

MARVEL 

dry  alluviiun 

DUG, 

TUC 

4.26 

2.2 

01/26/68 

CABRIOLET 

dry  rhyolite 

DUG, 

TUC 

4.30 

2.3 

03/12/68 

BUGGY  I 

dry  basalt 

BKS, 

DUG, 

TUC 

4.61 

5.4 

03/14/68 

POMMARD 

dry  alluvium 

DUG, 

tik: 

3.91 

1.4 

04/26/68 

BOXCAR 

rhyolite 

BKS 

6.62 

6.2 

1300 

12/08/68 

SCHOONER 

dry  tuff 

BKS, 

DUG, 

TUC 

5.31 

4.8 

35 

12/19/68 

BENHAM 

wet  tuff 

BKS 

6.65 

6.3 

1100 

10/29/69 

CALABASH 

tuff 

BKS, 

mJG 

5.88 

5.6 

110 

03/26/70 

HANDLEY 

wet  tuff 

BKS, 

DUG 

6.58 

6.4 

>1000 

05/26/70 

FLASK 

tuff 

BKS, 

DUG 

5.86 

5.5 

105 

10/17/70 

CARPETBAG 

wet  tuff 

BKS, 

DUG 

6.14 

5.8 

220 

O7/O8/71 

MINIATA 

tuff 

BKS, 

DUG, 

TUC 

5.93 

5.5 

80 

04/26/72 

STARWORT 

wet  tuff 

BKS, 

DUG, 

TUC 

5.85 

5.6 

85 

09/26/72 

DELPHINIUM 

dry  alluvium 

BKS, 

TOC 

4.87 

4.4 

15 

From  Sprlnfier  and  Kinnaman  [l971,  1975 J 


TABLE  4.  Comparison  of  seismic  moment  of  tectonic  release,  and  m^(Lg) 
deviations  from  mean  curve  (eq.  5)  for  NTS  e5{plosions  ^ 


Name 

Date 

M  X  10^^ 

0 

( newton- meters ) 

BENHAM 

12/19/68 

5.6 

+0.05 

BOXCAR 

04/ 26/ 68 

1.4 

-0.03 

GREELEY 

12/20/66 

3.1 

-0.03 

HALF  BEAK 

06/30/66 

1.0 

+0.05 

HANDLEY 

03/26/70 

2.4 

+0.01 

SCOTCH 

05/23/67 

0.3 

+0.08 

values  from  Wallace  ^1985] 


TABLE  5.  3-nd  m^(Lg)  values  of  NTS  exjjlosions  of  unannounced  yields 


Date 


Name 


Rock  Type 


m^(P-ISC)  m^(Lg) 


I2/05M 

04/21/65 
04/25/66 
12/13/66 
01/19/67 
02/08/6? 
02/23/67 
02/23/67 
05/10/67 
06/26/67 
06/29/67 
07/27/67 
08/1 8/67 
08/31/67 
09/27/67 
10/18/67 

II/O8/6? 

02/21/67 

02/29/68 

03/22/68 

08/29/68 

09/06/68 

09/17/68 

09/24/68 

11/20/68 

01/15/69 

01/30/69 

03/20/69 

03/21/69 

05/07/69 

05/27/69 

06/12/69 

07/16/69 

09/16/69 

10/08/69 

11/21/69 

12/05/69 

12/17/69 

12/17/69 

12/18/69 

02/04/70 

02/11/70 

02/25/70 

02/26/70 

03/23/70 

04/21/70 

04/21/70 


CREPE 
GUMDROP 
PIN  sraiPE 
NEW  POINT 
NASH 
WARD 

PERSIMMON 

AGILE 

MICKEY 

MIDI  MIST 

UMBER 

STANLEY 

BORDEAUX 

DOOR  MIST 

ZAZA 

LANPHER 

COBBIER 

KNOX 

DORSAL  FIN 

STINGER 

SLED 

NOGGIN 

STODDARD 

HUDSON  SEAL 

MING  VASE 

WINESKIN 

VISE 

BARSAC 

COFFER 

PURSE 

TORRIDO 

TAPPER 

ILDRIM 

JORUM 

PIPKIN 

PICCALILLI 

DIESEL  TRAIN 

LOVAGE 

GRAPE  A 

TERRINE 

GRAPE  B 

DIANA  MIST 

CUMARIN 

YANNIGAN 

SHAPER 

SNUBBER 

CAN 


tuff 

tuff 

tuff 

alluviiom 

dolomite 

alluvium 

alluvium 

tuff 

tuff 

tuff 

alluvium 

tuff 

alluvium 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

alluvium 

alluvium 

alluvium 

tuff 

tuff 

alluvium 

tuff 

tuff 

tuff/ rhyolite 

tuff 

tuff 

alluvium 

tuff 

tuff 

tuff 

tuff 

tuff 

alluvium 

tuff 

tuff 

tuff 


4.8 

5.14 

5.0 

5.06 

4.5 

5.10 

4.6 

4.96 

5.3 

5.51 

4.6 

4.77 

4.4 

4.25 

5.6 

5.91 

4.9 

5.26 

5-1 

5.20 

4.6 

4.82 

5.0 

5.19 

4.6 

5.01 

5.0 

5.05 

5-7 

6.00 

5.7 

6.09 

5-1 

5.10 

5.8 

6.01 

5.0 

5.14 

5.6 

6.00 

5-9 

6.12 

5.5 

6.08 

5.1 

5.24 

5.0 

5.22 

4.9 

5.08 

5.3 

5.58 

4.9 

5.23 

4.4 

4.73 

4.9 

5.15 

5.5 

5.98 

5.0 

5.44 

4.5 

4.85 

4.6 

5.28 

6.1 

6.52 

5.6 

6.06 

5.0 

5.26 

4.9 

5.23 

4.7 

4.96 

5.4 

5.79 

5.2 

5.60 

5.6 

6.0c 

4.7 

5.02 

5.2 

5.3^ 

5.3 

5.58' 

5.5 

5.69 

4.4 

4.91 

4.6 

5.13 

TABLE  5.  (continued) 


05/01/70 

05/05/70 

05/15/70 

05/21/70 

11/05/70 

12/18/70 

06/24/71 

09/21/72 

03/08/73 

04/25/73 

06/05/73 

06/06/73 

06/28/73 

10/12/73 

01/03/76 

03/14/76 


HOD 

MINT  LEAF 

CORNICE 

MORRONES 

ABEYTAS 

BANEBERRY 

HAREBELL 

OSGURO 

MIERA 

ANGUS 

DIDO  QUEEN 

ALMENDRO 

PORTUUCA 

HUSKY  ACE 

MUENSTER 

COLBY 


Rock  Type 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

tuff 

alluvium 

tuff 

tuff/rhyolite 

alluvium 

tuff 


m  (P-ISC)  nv^(Lg) 


5.29 

5.56 

5.44 

5.34 

5.21 

5-49 

5.97 

5.75 

5.09 

5.12 

6.31 

5.32 
4.83 
6.53 
6.53 


JANUARY  19,  1967,  I6"45'"00 


ALMENORO,  JUNE 


-»**'-*  »”►  k**  k"*  VNV  'k  •  L  •  L'^  .  • 


HANDCAR 


O  SHOAL 


DISCUS  THROWER 


9REELEY 


INTER-STATION  SURFACE  WAVE  ANALYSIS  BY 


FREQUENCY-DOMAIN  WIENER  DECONVOLUTION  AND  MODAL  ISOLATION. 

By 

Jay  Horng-Jye  Hwang  and  B.  J.  Mitchell 

ABSTRACT 

A  new  technique  which  combines  frequency-domain  Wiener  filtering 
and  modal  isolation  is  developed  for  determining  interstation  phase  veloci¬ 
ties,  group  velocities,  and  attenuation  coefficients  of  seismic  surface  waves. 
Frequency-domain  Wiener  filtering  is  more  effective  than  time-domain 
Wiener  filtering  for  the  determination  because  it  uses  a  smaller  window 
lag  which  produces  a  smoother  interstation  Green’s  function.  This  leads 
to  greater  accuracy  and  stability  when  noise-contaminated  data  are 
analyzed.  We  optimize  Wiener  filtering  in  the  frequency  domain  by 
applying  two  trapezoidal  windows  of  different  lags  to  the  cross-correlation 
function  between  two  stations  and  to  the  autocorrelation  function  of  the 
first  station,  respectively.  The  windowed  correlation  functions  are  then 
transformed  to  the  frequency  domain.  The  interstation  Green’s  function 
in  this  technique  is  the  ratio  of  the  smoothed  cross-spectrum  to  the 
smoothed  autospectrum  of  the  first  station.  Frequency-domain  Wiener 
filtering  is  equivalent  to  time-domain  Wiener  filtering  when  the  same  rec¬ 
tangular  window  is  applied  to  the  both  correlation  functions. 

Wiener  filtering,  however,  cannot  efficiently  remove  higher  mode 
interference  when  the  higher  modes  are  superimposed  on  the  fundamental 
mode  in  the  correlation  functions.  To  more  thoroughly  eliminate  the 


effects  of  such  interference,  phase- matched  filtering  or  time-variable 
filtering  can  be  employed  to  isolate  one  particular  mode  at  each  of  two 
stations.  Frequency-domain  Wiener  deconvolution  is  then  applied  to  cal¬ 
culate  the  Green’s  function.  The  interstation  group  velocities  can  be 
obtained  by  applying  the  multiple  filter  technique  to  the  Green’s  function, 
and  can  be  refined  by  phase-matched  filtering.  The  amplitude  and  phase 
spectra  of  the  Green’s  function  are  used  to  calculate  attenuation 
coefficients  and  phase  velocities,  respectively,  for  the  interstation  medium. 

This  new  technique  is  compared  with  other  methods  by  applying 
them  to  both  noise-contaminated  synthetic  seismograms  and  real  data. 

The  proposed  technique  is  found  to  be  superior,  particularly  in  period 
ranges  where  the  signal-to-noise  is  low. 

INTRODUCTION 

In  surface  wave  studies,  seismologists  attempt  to  determine  group  arrival  time, 
phase  angle,  and  amplitude  as  a  function  of  period.  These  quantities  are  important  for 
studies  of  the  crust  and  upper  mantle  structure,  earthquake  source  mechanisms,  and 
anelastic  properties  of  the  earth.  Determinations  of  these  quantities  have  commonly  been 
done  using  either  single-station  or  two-station  methods. 

Single-station  methods  assume  either  that  the  initial  phase  of  the  earthquake  source 
is  known  or  that  uncertainty  arising  from  lack  of  knowledge  of  the  initial  phase  is  small 
enough  to  be  ignored.  If  that  is  the  case,  methods  such  as  the  moving-window  analysis 
(Landisman  etal,  1969)  or  multiple-filter  technique  (Dziewonski  etal,  1969)  can  be  used 
for  group  velocity  determinations  and  observed  phases  can  be  used  to  determine  phase 
velocity.  Knopoff  and  Schwab  (1968)  have  shown,  however,  that  corrections  to  these 
measurements  are  required  whenever  the  seismic  source  is  not  oriented  in  either  a  purely 
horizontal  or  a  purely  vertical  direction.  These  corrections  can  be  substantial,  especially 


at  short  distances  from  the  source. 

A  single-station  method  has  also  been  developed  for  inferring  Q  structure  of  the 
crust,  from  the  surface  wave  amplitudes  (Cheng  and  Mitchell,  1981).  This  method,  like¬ 
wise,  requires  a  knowledge  of  the  earthquake  source  mechanism.  All  single-station 
methods,  therefore,  require  prior  knowledge  of  the  earthquake  focal  mechanism  or  an 
assumption  that  its  effects  are  small  enough  to  be  ignored  when  they  are  used  to  deter¬ 
mine  surface  wave  velocities  or  amplitudes. 

Two-station  methods  avoid  the  necessity  of  knowing  the  earthquake  focal  mechan¬ 
ism,  but  they  require  two  seismic  stations  which  lie  on  a  common  great  circle  path  with 
the  earthquake.  Such  pairs  of  stations  can  be  used  to  measure  interstation  attenuation 
coefficients  by  taking  spectral  ratios,  and  to  measure  interstation  phase  velocity  by  the 
Fourier  phase  method  (Sato,  1955,1956).  Interstation  group  velocity  can  be  calculated 
by  multiply  filtering  each  seismogram,  and  then  dividing  the  group  delay  into  the  station 
separation. 

Landisman  et  al.  (1969)  suggested  that  by  windowing  a  cross  correlogram  one  can 
reduce  noise  and  stabilize  the  measured  phase  velocity.  They  also  noted  that  the  cross¬ 
correlation  function  can  approximate  the  interstation  impulse  response,  and  that  moving 
window  analysis  or  multiple-filter  technique  (MFT)  can  be  used  to  obtain  the  intersta¬ 
tion  group  velocities.  Nakanishi  (1979)  presented  a  Wiener  filtering  technique  to  measure 
phase  velocity  and  Q  using  Rayleigh  waves  on  a  great  circle  path  recorded  at  single  sta¬ 
tion.  More  recently,  Taylor  and  Toksdz  (1982)  used  a  time-domain  Wiener  filtering  tech¬ 
nique  to  obtain  interstation  phase  and  group  velocities  and  attenuation  coefficients. 
Time-variable  filtering  was  introduced  by  Pilant  and  Knopoff  (1964)  to  smooth  surface 
wave  trains  contaminated  by  beats.  Landisman  el  al.  (1969)  used  this  filter  to  isolate  a 
single  surface  wave  mode  before  calculating  the  interstation  phase  velocity.  More 
recently,  phase-matched  filtering  has  been  used  to  extract  the  primary  wave  train  from 
the  multipathing  event  (Herrin  and  Goforth,  1977). 


In  this  paper,  we  propose  a  new  technique  which  combines  signal  isolation  and 
frequency-domain  Wiener  filtering  to  calculate  the  interstation  Green’s  function  for  sur¬ 
face  waves.  The  process  of  frequency-domain  Wiener  filtering  will  first  be  presented  and 
the  results  of  its  applications  will  be  compared  with  that  of  time-domain  Wiener  filtering 
and  the  spectral  ratio  technique.  Signal  isolation  using  phase-matched  or  time-variable 
filtering  will  be  discussed  in  the  following  section.  Synthetic  seismograms  will  be  used  to 
test  frequency-domain  Wiener  filtering  against  other  techniques,  and  then  used  to  test 
the  effectiveness  of  isolation  by  phase-matched  and  time-variable  filtering.  The  proposed 
technique  will  also  be  applied  to  the  real  data. 

INTERSTATION  GREEN’S  FUNCTION 

The  interstation  Green’s  function  corresponding  to  the  medium  impulse  response 
can  be  estimated  if  there  are  two  seismograms  positioned  along  the  same  great  circle 
path  from  a  source.  The  seismogram  at  station  I,  the  nearer  station  to  the  source,  can 
be  treated  as  the  input  to  the  interstation  crustal  medium.  The  input  propagates 
through  the  interstation  crust  and  produces  the  output  recorded  at  station  2.  Let  x, 
represent  the  input  to  the  system  h„  and  y^  be  the  output  where 

Xt  =  ,  X,  ,  .  ,  x„  )  (1) 

fit  =  (fio  >  fit  >  fi2 .  .  fim  ) 

yi  =  (yo ,  yi ,  yj , . ,  y,+m) 

The  relationship  among  them  can  be  expressed  by 

r*in 

yt  =  S  fir  Xt-r  for  t=0,l,2, .  ,n-)-m  (2) 

f^O 

in  the  time  domain  and 

X(f)  H(f)  =  Y(f)  (3) 

in  the  frequency  domain  where  X(f),  H(f),  and  Y(f)  are  Fourier  spectra  of  total  lengths  of 
X,  ,yi  and  h,  .  The  goal  is  to  determine  the  Green’s  function  or  filter  h,.  There  are  various 
deconvolution  schemes  can  be  used  to  solve  for  it. 


SPECTRAL  RATIO 


From  equation  (3),  the  input  can  be  easily  deconvolved  by  dividing  (3)  by  the  input 
and  computing  the  Green’s  function 


H(f)  = 


m  =  -LYIO  I 

X(f)  I  X(f)  1 


(4) 


where  and  4>y  are  phase  spectra  of  X(f)  and  Y(f),  respectively.  This  simple  deconvolu¬ 
tion  could  be  unstable,  particularly  in  the  presence  of  spectral  holes.  When  computing 
X(f)  and  Y(f),  the  total  lengths  of  Xi  and  y,  are  used.  This  technique  thus  considers  that 
all  of  the  points  are  signal.  Since  this  method  cannot  distinguish  signal  of  individual 
modes  from  other  disturbances,  the  deconvolution  results  could  be  inaccurate  because  of 
the  random  noise,  multipathing,  or  interference  by  other  modes  which  are  usually 
present  in  real  data. 


TIME-DOMAIN  WIENER  FILTERING 


A  least-squares  or  Wiener  deconvolution  (Wiener  1949;  Treitel  and  Robinson,  1966; 
Peacock  and  Treitel,  1969)  has  been  employed  by  Taylor  and  Toksoz  (1982)  on  surface 
waves  to  obtain  the  impulse  response  of  the  interstation  medium  (Taylor  and  Toksoz, 
1982).  Using  this  scheme,  the  problem  is  to  determine  the  h,  such  that  the  actual  output 
signal  Yi  is  the  least-squares  approximation  to  the  desired  output  z, 

—  (*oi*l>*2> . *o+iii)  (5) 

The  error  between  desired  and  actual  output  is 

et  =  yt  -  z,  for  t=0,l,2, . n-(-m  (6) 


From  equations  (2)  and  (6),  the  equation  to  be  solved  in  matrix  form  is 
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In  matrix  shorthand, 


Xh=y+e  (7b) 

From  the  theory  of  least  squares,  the  normal  equations  in  matrix  form  are 

X'rXh  =  X'^y  (8) 

where  X^  is  the  transpose  of  X.  It  turns  out  that  terms  X^X  give  the  autocorrelation  of 

input  X,.  X  ^X  can  be  written 

=  E  ’^t+r  fo*"  *i=0.1,2, . m  (9) 

r=0 

and  X  give  the  cross-correlation  function  of  input  and  output 

gi  =  E  Xryt+r  for  t=0,l,2, . m  (10) 

r*»0 

Thus  the  m-length  Wiener  filter  can  be  obtaii.vd  from  the  solution  of  the  normal  equa¬ 
tions  of  the  form 


Since  the  autocorrelation  matrix  in  equation  (11)  is  in  Toeplitz  form,  the  Levinson  recur¬ 
sion  (Wiener  1949;  Treitel  and  Robinson,  1966)  can  be  used  to  solve  this  equation  with 
high  efficiency  and  minimum  of  computer  storage.  However,  numerical  instability  might 
occur  for  large  m.  To  stablize  the  deconvolution,  Taylor  and  Toksds  (1982)  added  a  small 
constant  6'^  to  the  autocorrelation  at  zero  lag,  i.e.,  to  the  diagonal  elements  in  equation 
(11).  The  equation  was  thus  solved  with  a  damped  least-squares  technique.  Since  damp¬ 
ing  has  altered  the  filter  estimate,  the  recovery  of  true  solution  can  be  performed  in  fre¬ 
quency  domain  by  scaling  the  altered  filter  estimate 

H(r)  =  H'  (f)  (12) 

where  H'  (f)  is  the  filter  estimated  by  a  damped  least-squares  technique. 

FREQUENCY-DOMAIN  WIENER  FILTERING 
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In  lime-domain  Wiener  filtering,  the  normal  equations  are  formed  in  terms  of  the 
autocorrelation  and  cross- correlation  functions,  and  the  deconvolution  scheme  is  carried 
out  in  the  time  domain.  The  deconvolution  can  be  considerably  simplified  if  the  normal 
equations  in  equation  (11)  are  transformed  to  the  frequency  domain  (Jenkins  and  Watts, 
1968).  Note  that  there  are  only  m  points  taken  in  both  correlation  functions  in  time- 
domain  Wiener  deconvolution.  It  is  equivalent  to  windowing  both  complete  correlation 
(i.e.  n-fm  points)  functions  with  one  rectangular  window  of  lag  m.  Since  windowing  in 
the  time  domain  results  in  smoothing  in  the  frequency  domain,  the  spectra  from 
transforming  m-length  correlation  functions  are  smoothed  spectra  instead  of  the  original 
spectra  which  are  the  Fourier  transform  of  n-|-m-length  correlation  functions.  Thus 
Fourier  transforming  equation  (11)  gives 

G(f)  =  H(f)  R(f)  (13) 

where  G(f),  H(f),  and  R(f)  are  Fourier  spectra  of  g„  h,,  and  Tj  of  length  m,  respectively. 

The  deconvolution  becomes  very  straightforward 


H(f)  =  SSL 
R(f) 

Therefore  the  Green’s  function  in  the  frequency  domain  is  actually  just  the  ratio  of  the 
smoothed  cross-spectrum  to  the  smoothed  autospectrum.  It  avoids  dealing  with  the  ins¬ 
tability  and  necessary  correction  in  time-domain  Wiener  filtering. 

At  this  point,  it  is  apparent  that  the  advantage  of  Wiener  filtering  over  the  spectral 
ratio  technique  is  windowing  which  cuts  off  the  noise  outside  the  window.  However,  if 
we  note  that  the  autocorrelation  function  usually  has  a  shorter  duration  than  the  cross 
correlation  function  because  the  autocorrelation  has  zero  phase,  we  can  improve  the 
Wiener  filtering  by  applying  a  shorter  window  to  the  autocorrelation  function  to  achieve 
a  larger  reduction  of  noise.  In  addition,  since  the  rectangular  window  has  a  big  side  lobe 
which  results  in  leakage  when  convolution  is  performed,  we  instead  use  a  trapezoidal 
window  which  is  composed  of  rectangular  and  Parzen  windows  (Jenkins  and  Watts, 
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where  and  tp  indicate  the  lags  of  the  rectangular  and  Parzen  windows,  respectively. 
The  Parzen  window  is  used  to  reduce  the  sharpness  of  the  discontinuity  of  the  rectangu¬ 
lar  window.  Thus  the  proposed  frequency-domain  Wiener  filtering  is  summarized  as  fol¬ 
lows:  (1)  The  cross  correlation  and  autocorrelation  functions  of  the  total  seismograms 
duration  are  calculated,  and  are  windowed  by  trapezoidal  windows  of  different  lags.  The 
lags  are  visually  optimized  to  achieve  the  largest  possible  noise  reduction,  (2)  the  win¬ 
dowed  correlation  functions  are  Fourier  transformed  into  the  frequency  domain  to  obtain 
the  smoothed  spectra,  and  (3)  The  Green’s  function  is  simply  the  ratio  of  the  smoothed 
cross  spectrum  to  the  smoothed  autospectrum. 


TESTS 

Synthetic  Rayleigh  seismograms  are  generated  at  distances  of  1000  km  and  2000 
km  from  the  source.  Tests  are  run  for  both  a  noiseless  case  and  for  seismograms  contam¬ 
inated  by  additive  Gaussian  random  noise.  In  both  cases,  phase  and  group  velocities,  and 
attenuation  coefficients  are  calculated  from  Green’s  functions  by  the  three  techniques 
mentioned  above.  The  interstation  group  velocities  are  calculated  by  applying  the  MFT 
to  the  Green’s  functions,  and  the  interstation  phase  velocities  are  calculated  from  the 
phase  spectra  of  transfer  functions  using  the  formula 
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where  Ad  is  the  interstation  distance,  to  is  the  first  time  point  of  the  Green’s  function 
and  <^{f)  is  the  phase  of  the  Green’s  function  in  cycles.  Attenuation  coefficients  of  the 
interstation  medium  were  calculated  using 
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where  A,  and  Aa  are  the  epicentral  distances  in  degrees  of  station  1  and  station  2. 
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The  model  used  to  generate  synthetic  Rayleigh  seismograms  is  listed  in  Table  1. 
The  seismograms  are  shown  in  Figure  la.  Three  different  deconvolution  schemes  are 
applied,  and  the  results  are  shown  in  Figure  lb.  Not  surprisingly,  the  phase  and  group 
velocities,  and  the  attenuation  coefficients  by  the  three  techniques  all  agree  exactly  with 
the  theoretical  values.  Since  there  is  no  noise  contamination,  the  windowing  in  the 
Wiener  filtering  process  does  not  produce  any  improvement. 

Gaussian  random  noise  is  then  generated,  which  has  zero  mean  and  for  which  the 
standard  deviation  is  in  fraction  of  mean  absolute  amplitude  of  the  synthetic  seismo¬ 
gram.  In  the  first  test,  random  noise  with  an  amplitude  which  is  30%  of  mean  absolute 
amplitude  of  the  signal  is  added  to  the  synthetic  seismograms  in  Figure  la.  The  noise- 
contaminated  seismograms  and  spectra  are  shown  in  Figure  2a.  Note  that  the  S/N 
varies  with  period.  It  is  larger  in  the  intermediate  period  range,  and  smaller  at  longer 
and  shorter  periods. 

An  interactive  program  was  written  to  view  and  shift  the  correlation  functions  so 
that  the  lags  of  the  window  can  be  optimized  to  achieve  as  large  a  noise  reduction  as 
possible.  Figure  2b  shows  the  subtle  difference  between  the  time-domain  Wiener  filtering 
and  the  frequency-domain  Wiener  filtering.  The  lag  of  the  window  applied  to  the  auto¬ 
correlation  function  in  frequency-domain  Wiener  filtering  is  smaller  than  the  lag  used  in 
time-domain  Wiener  filtering,  while  the  lags  of  the  window  applied  the  cross-correlation 
function  are  the  same  in  both  techniques.  The  results  in  the  tests  using  frequency- 
domain  and  time-domain  Wiener  filtering  (not  shown  here)  are  exactly  the  same,  as  long 
as  the  rectangular  windows  applied  to  both  correlation  functions  have  the  same  lag. 

The  Green’s  functions  determined  by  different  techniques  are  shown  in  Figure  2c 
where  it  is  seen  that  the  Green’s  function  obtained  by  frequency-domain  Wiener  filtering 
is  smoother  and  cleaner  than  that  by  time-domain  Wiener  filtering.  The  Green’s  function 
by  the  spectral  ratio  technique  is  the  poorest  of  the  three,  because  it  cannot  reduce  level 
of  the  noise.  Figure  2d  shows  group  and  phase  velocities  and  attenuation  coefficients.  As 


in  noiseless  case,  the  phase  and  group  velocities  by  different  schemes  are  equally  correct. 
Small  errors  due  to  random  noise  contamination  appear  in  group  velocity,  but  not  in 
phase  velocity,  indicating  that  the  phase  velocity  is  less  sensitive  than  group  velocity  to 
the  noise.  However,  the  attenuation  coefficients  calculated  from  the  three  techniques  are 
quite  scattered,  especially  at  the  long  periods  The  spectral  amplitudes  of  the  Green’s 
function  at  long  periods  is  very  close  to  unity,  as  shown  in  Figure  lb.  When  the  loga¬ 
rithmic  value  is  taken  as  in  equation  (17),  a  slight  contamination  on  the  Green’s  function 
will  result  in  a  large  error  in  the  attenuation  coefficient.  For  example,  the  logarithmic 
value  of  0.965  for  H(f)  leads  to  a  attenuation  coefficient  value  which  is  1.4  times  that  for 
which  the  logarithmic  value  of  H(f)  is  0.975  This  means  that  a  1%  of  error  in  the 
Green’s  function  can  lead  to  40%  of  error  in  attenuation  coefficient.  For  this  reason,  the 
attenuation  coefficients  at  long  periods  are  very  sensitive  to  the  noise  contamination. 
The  results  of  Figure  2d  show  that  the  attenuation  coefficients  determined  by  the 
frequency-domain  Wiener  filtering  are  more  accurate  and  more  stable  than  those  by 
other  two  techniques.  This  occurs  because  frequency-domain  Wiener  filtering  is  more 
capable  of  reducing  noise  contamination,  and  provides  a  smoother  and  more  accurate 
estimate  of  the  interstation  response. 

In  real  surface  wave  data  interference  phenomena  often  occur,  due  to  overlap 
among  modes  or  to  multipathing  effects.  To  examine  the  windowing  effect  in  frequency- 
domain  Wiener  filtering  for  these  cases,  synthetic  seismograms  of  Rayleigh  waves  includ- 
iig  fundamental  and  first  higher  mode  were  generated  from  the  same  model  in  Table  1 
for  distances  of  1000  km  and  2000  km.  Gaussian  random  noise  with  zero  mean  and  stan¬ 
dard  deviation  of  20%  of  mean  absolute  amplitude  are  added  to  the  seismograms.  The 
seismograms  and  spectra  are  shown  in  Figure  3a.  Spectral  ratio  and  frequency-domain 
techniques  are  applied  to  analyze  the  data.  The  Green’s  functions  are  shown  in  Figure 
3b.  It  can  be  seen  that  windowing  did  not  remove  the  interference  effect  caused  by  the 
higher  mode.  The  reason  that  windowing  did  not  work  becomes  apparent  from  the 


correlation  functions  in  Figure  3c.  It  is  mainly  because  the  first  higher  mode  rides  on  the 
fundamental  mode  in  the  cross-correlation  function,  thus  cannot  be  removed  by  window¬ 
ing.  Therefore  it  is  desirable  to  remove  the  first  higher  mode,  i.e.,  to  extract  the  funda¬ 
mental  mode  before  Green’s  function  is  calculated.  In  next  section,  the  time-variable  and 
phase-matched  filters  are  applied  for  this  purpose.  The  MFT  is  also  used  for  comparison. 

TIME- VARIABLE  FILTERING 

Time-variable  filtering  allows  us  to  use  a  given  dispersion  curve  to  extract  a  wavetrain 
from  a  complex  signal.  Such  filtering  is  necessary  when  the  phase  or  amplitude  of  the 
extracted  wavetrain  has  to  be  analyzed.  Given  a  group  velocity  curve,  « (f),  correspond¬ 
ing  to  signal,  /  (t),  this  dispersion  data  may  be  written  as  the  group  delay  at  frequency, 

f. 
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where  d  is  the  epicentral  distance  in  kilometers.  Let  F(f)  be  the  Fourier  spectrum  of 
/  (t),  i.e., 
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/  (t)  = /F(f)e‘2’'f‘df  (19) 
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This  equation  means  that  the  Fourier  synthesis,  i.e.,  the  interference  among  all  the  com¬ 
ponent  harmonic  waves  results  in  /  (t).  The  dispersion  curve  implies  that  the  harmonic 
waves  interfere  constructively  in  the  vicinity  of  the  curve  tg(f).  Thus  isolation  can  be 
done  by  confining  the  interferences  to  the  vicinity  of  the  curve  tg(f).  Filtering  requires 
specification  of  the  group  arrival  times  corresponding  to  the  beginning  and  end  of  the 
transmission  window  for  each  of  the  component  waves.  Outside  the  transmission 
domain,  we  can  neglect  the  contribution  of  Fourier  components  to  the  signal  we  intend 
to  isolate.  Therefore  we  can  resynthesize  signal  which  pertains  to  the  dispersion  curve, 
tt(f),  by  summing  the  Fourier  components  truncated  at  the  upper  and  lower  limits  of  the 
transmission  domain.  The  isolated  signal  is 


00 

/  (t)  =  /  F(f)  w  (t.f)  e‘2^ '  ‘  df  (20) 
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where  w(f,t)  is  a  window  which  is  zero  outside  the  transmission  domain.  In  equation 
(20),  multiplication  of  a  spectral  line  F(fo)  by  a  time  window  w(t,fo)  results  in  spreading  of 
this  line  into  a  spectrum  of  finite  width,  i.e. 
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w  (  fo,f-fo)  =  /  w  (  t  ,  fo  )  dt  (21) 
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The  spectrum  of  the  filtered  seismogram  can  be  expressed 

00 

F(fo)=  /F(f)W(U„-f)dr  (22) 
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which  simply  represents  convolution  of  spectra  F(f)  and  W(f).  Therefore,  the  spectrum 
of  the  isolated  spectrum  is  also  a  smoothed  spectrum  of  the  original  spectrum.  It  is 
apparent  from  equation  (22)  that  the  performance  of  the  time-variable  filtering  relies 
totally  upon  W(f).  The  choice  of  the  window  is  thus  very  important. 

Jenkins  and  Watts  (1968)  pointed  out  that  smoothing  reduces  the  variance  ,  i.e., 
increases  the  smoothness  but  increases  the  bias  of  the  estimator.  A  sensible  procedure  is 
to  compromise  by  making  the  sum  of  variance  and  bias  as  small  as  possible.  Various 
types  of  windows  are  adopted  in  time-variable  filtering  to  check  the  smoothness  and 
biases  of  the  isolated  spectra.  We  found  that  smoothing  effect  of  the  windows  varied 
from  Gaussian,  Parzen,  Tukey,  and  cosine  windows,  in  that  order,  with  Gaussian  win¬ 
dow  exhibiting  the  greatest  smoothing.  The  biases,  however,  are  in  reverse  order.  Our 
numerical  experiments  show  that  the  Tukey  window  is  best  for  signals  with  a 
moderately-peaked  spectrum.  The  fundamental  mode  often  exhibits  such  a  spectrum. 
The  Tukey  window  is  given  by: 


Wt ( t )  = 


— (  14-cos - 1 


t-tj  I  <  t. 


>  W 


(23a) 


where  t^  indicates  the  lag  of  the  window.  For  a  highly-peaked  spectrum  such  as  that 
often  pertaining  to  higher  modes,  the  cosine  window  is  adopted  because  it  has  the  smal- 
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lest  biasing  effect.  For  such  a  highly-peaked  spectrum  smoothing  will  lead  to  very 
significant  biases  to  which  the  estimates  of  attenuation  coefficients  are  very  sensitive. 
The  cosine  window  is  given  by: 

Wc  ( t )  =  ,  I  t-tg  I  <  t* 

1  I  ^  ^ 

After  the  type  of  the  window  is  decided,  the  remaining  question  left  is  the  choice  of  the 
window’s  lag.  Landisman  et  al.  (1969)  used  for  t*: 

t«  =  T  (  «  +  /3  -^  )  (24) 

dti 

where  T  is  a  period,  a,  are  empirical  parameters.  We  find  that  the  term  leads  to 

fluctuating  behavior  of  the  spectrum,  so  0  is  taken  to  be  zero  in  the  present  study.  Cara 
(1973)  optimized  the  lag  of  the  window  by  solving  the  equation  relating  the  duration  of 
the  filtered  signal  to  the  bandwidth  of  the  Gaussian  filter  and  the  dispersiveness  of  the 
wavetrain.  We  have  adopted  an  empirical  approach,  suggested  by  Jenkins  (1961),  which 
uses  smaller  lags  initially  and  then  progressively  larger  lags.  By  allowing  the  lag  to 
become  larger,  more  significant  detail  of  the  spectrum  can  be  explored,  because  the 
bandwidth  of  W  is  narrower.  In  this  application,  the  process  of  increasing  the  lag  contin¬ 
ues  until  the  spectral  biases  do  not  result  in  significant  systematic  errors  in  7.  This 
empirical  procedure  was  later  referred  to  as  window  closing  (Jenkins  and  Watts,  1968). 

The  MFT  is  applied  to  provide  the  group  velocity  for  the  time-variable  filtering  to 
calculate  the  group  delay.  Following  the  window  closing  procedure,  we  find  that  the 
Tukey  window  with  a  time  lag  with  four  or  five  times  the  period  is  large  enough  to  give 
a  correct  spectral  estimate.  To  obtain  the  Green’s  function,  the  fundmental  mode  of 
each  of  the  seismograms  in  Figure  3a  is  extracted  first,  then  frequency-domain  Wiener 
filtering  is  applied.  The  seismogram  and  spectra  at  2000  km  in  Figure  3a  are  shown 
before  and  after  isolation  in  Figure  4a.  It  can  be  seen  that  most  of  the  higher  mode 
interference  and  noise  contamination  have  been  removed.  The  group  and  phase  veloci- 


ties,  and  attenuation  coefficients  in  Figure  4b  agree  well  with  theoretical  values.  It 
should  be  noted,  however,  that  the  spectral  ratio  technique  works  almost  as  well  as 
frequency-domain  Wiener  filtering  after  isolation,  since  there  is  little  noise  left. 

The  amplitudes  of  individual  modes  at  each  station  can  also  often  be  discriminated 
if  they  are  isolated  on  MFT  plots.  In  another  approach,  the  amplitude  spectra  of  higher 
modes  with  nearly  identical  group  velocities,  can  be  treated  as  a  superposition  of  modes 
and  compared  with  theoretical  components  of  higher  modes.  This  approach,  called  the 
multi-mode  method,  was  developed  by  Cheng  and  Mitchell  (1981).  It  has  been  used  to 
study  regional  variations  of  crustal  anelasticity  along  relatively  short  paths.  In  the  exam¬ 
ple  which  follow,  We  will  consider  only  a  single  higher  mode,  but  the  method  will  be 
equally  applicable  to  any  number  of  isolated  modes  or  to  multi-mode  analysis  if  the 
higher  modes  can  not  be  separated.  The  attenuation  coefficients  of  the  different  modes 
can  be  calculated  by  taking  the  spectral  ratio  of  the  amplitudes  from  MFT  in  the  same 
way  that  fundamental  mode  are  obtained.  -7  values  calculated  from  the  spectral  ampli¬ 
tudes  by  the  MFT  are  also  presented  in  Figure  4c.  The  7  values  determined  from  MFT 
amplitudes  are  not  as  smooth  or  as  accurate  as  those  in  Figure  4b  determined  by  the 
proposed  technique.  The  systematic  errors  in  7  in  the  period  range  from  20  sec  to  40  sec 
are  caused  by  the  large  spectral  amplitude  biases  produced  by  Gaussian  windowing.  In 
isolating  the  higher  mode,  a  cosine  window  with  a  lag  of  five  times  the  period  is  used  in 
the  time-variable  filter.  The  spectral  biases  of  higher  mode  obtained  by  MFT  (Figure 
4d)  are  very  significant,  although  those  of  fundamental  mode  are  usually  less  significant. 
The  7  estimates  of  the  higher  mode  from  the  isolated  amplitudes  obtained  by  MFT  and 
time-variable  filtering  are  also  shown  in  Figure  4d.  At  a  period  of  9  sec,  the  error  in  the 
7  value  by  using  MFT  is  84%,  while  the  error  by  using  time-variable  filtering  is  8%.  At  a 
period  of  10  sec,  the  error  in  7  from  MFT  is  104%,  while  error  from  time-variable  filter¬ 
ing  is  22%.  We,  therefore,  that  the  technique  we  propose  in  this  paper  also  works  very 
well  for  higher  modes. 
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PHASE- MATCHED  FILTICRING 


A  phase-matched  filter  is  a  filter  in  which  the  Fourier  phase  of  the  filter  is  matched 
to  that  of  a  signal.  To  compromise  between  S/N  improvement  and  time  resolution,  a 
white  spectrum  is  adopted  for  the  filter  in  this  application.  The  phase  spectrum  needs  to 
be  obtained.  From  Papoulis  (1962)  the  group  delay  and  Fourier  phase  of  a  signal  are 
related  as  follows 

“n 

<()  (wo)  =  /  tg(w)dw  (25b) 

0 

where  u  is  the  angular  frequency.  tg(oj)  can  be  calculated  from  equation  (18)  in  which  the 
group  velocity  can  be  provided  with  the  MFT  analysis  on  the  wavetrain.  To  complete 
the  construction  of  the  filter,  the  Fourier  phase  of  the  filter  is  computed  from  equation 
(25b).  The  filter  is  then  correlated  with  the  wavetrain.  Since  the  filter  phases  are  close 
to  the  phases  of  the  signal,  the  phases  of  the  correlation  function  (so-called  pseudo¬ 
autocorrelation  function)  in  the  frequency  range  of  interest  are  small.  The  signal  on 
pseudo- autocorrelation  function  (PAF)  is  thus  spiked  so  that  we  can  apply  a  window  to 
reject  interfering  noise  outside  the  window,  i.e.,  isolate  the  signal  inside  the  window.  The 
phase  spectrum  of  the  windowed  PAF  is  then  used  to  correct  the  group  delay  of  the  trial 
filter  with  equation  (25a).  The  procedure  is  repeated  until  the  phase  spectra  of  the  filter 
and  the  signal  in  the  frequency  band  of  interest  are  identical.  The  phase-matched  filter 
can  therefore  be  employed  to  refine  the  group  velocity  and  to  extract  a  signal  from  a 
wavetrain.  As  in  Wiener  filtering,  the  effect  of  the  windowing  totally  relies  upon  the 
separation  of  noise  and  signal  in  the  correlation  function. 

We  applied  a  phase-matched  filter  to  extract  the  fundmental  mode  from  a  synthetic 
seismogram  at  a  distance  of  1000  km  in  Figure  3a.  The  trapezoidal  window  in  equation 
(15)  is  adopted.  The  PAF  and  the  spectral  before  and  after  isolation  are  shown  in  Fig¬ 
ure  5a.  We  can  see  that  higher  mode  interference  is  only  partially  removed,  and  the 


amplitudes  of  the  long  periods  are  not  well  determined.  This  is  because  the  separation 
between  the  fundamental  mode  and  higher  mode  on  the  PAF  is  not  large  enough.  The 
separation  is  larger  on  the  PAF  at  larger  distances;  in  such  cases  isolation  by  phase- 
matched  filtering  is  more  successful. 

Two  seismograms  and  their  spectra  shown  in  Figure  5b  are  generated  from  the 
same  model  at  distances  of  2000  km  and  4000  km.  Fundamental  and  first  higher  modes 
plus  random  noise  are  included.  Similarly  the  S/N  is  high  in  the  intermediate  period 
range,  while  low  in  the  long  and  short  period  ranges.  The  spectral  ratio  technique  is 
used  to  analyze  the  Green’s  function  after  the  isolation  by  phase-matched  filter.  The 
results  of  attenuation  coefficients  by  spectral  ratio  technique  with  and  without  modal 
isolations  by  phase-matched  filtering  are  shown  in  Figure  5c.  The  attenuation  coefficients 
obtained  by  using  time-variable  filter  isolation  are  also  shown  in  the  same  figure.  It  can 
be  seen  that  7  values  obtained  by  the  spectral  ratio  technique  without  modal  isolation 
are  quite  scattered  in  the  short  and  long  period  ranges  ,  since  the  S/N  is  low  in  these 
period  ranges.  The  qf  values  by  the  spectral  ratio  technique  with  isolation,  however,  show 
good  agreement  with  the  theoretical  values.  Phase  and  group  velocities  (not  shown)  also 
agree  very  well  with  the  theoretical  values.  Thus  the  proposed  technique  works  well  not 
only  in  high  S/N  situations,  but  also  does  when  S/N  is  low. 

APPLICATION 

The  proposed  technique  is  applied  to  a  surface  wave  path  from  an  event  occurring 
in  Alma-ata,  USSR,  on  18  October  1965  (10:21:43.8;  42.02N,  77.53E;  mi,=5.1  ;  depth=l 
km)  crossing  the  Indian  shield  between  WWSSN  stations  NDI  and  KOD  (6=13.31,  31.66, 
respectively).  Figure  6a  shows  the  vertical  component  seismogram  of  each  station  (after 
deconvolving  the  instrument)  before  and  after  isolation  for  fundamental  mode  by  time- 
variable  filtering.  The  interstation  Rayleigh  wave  phase  and  group  velocities  and 
attenuation  coefficients  are  shown  in  Figure  fib.  The  results  by  spectral  ratio  technique 


are  also  included  in  the  same  figure.  As  discussed  earli<'r,  the  amplitude  spectrum  of  the 
Green’s  function  at  long  periods  is  very  sensitive  to  noise  contamination.  This  figure 
shows  that  the  attenuation  coefficient  and  group  velocity  values  obtained  by  the  spectral 
ratio  technique  without  modal  isolation  are  very  scattered  at  long  periods.  Since  the  iso¬ 
lation  can  remove  the  effects  of  multipathing  and  higher  modes,  the  group  velocity  and 
attenuation  coefficients  by  the  proposed  technique  appear  to  be  smoother.  Phase  veloci¬ 
ties,  however,show  little  improvement,  since  they  are  less  sensitive  to  random  noise. 

CONCLUSION 

Time-domain  Wiener  filtering  is  found  to  be  nothing  more  than  the  ratio  of  the 
smoothed  cross-spectrum  to  the  smoothed  autospectrum  in  the  frequency  domain.  We 
have  optimized  Wiener  filtering  in  the  frequency  domain  by  applying  two  trapezoidal 
windows  of  different  lags  to  the  cross-correlation  and  autocorrelation  functions,  respec¬ 
tively.  Frequency-domain  Wiener  filtering  achieves  a  larger  reduction  of  noise  than  time- 
domain  Wiener  filtering.  Wiener  filtering,  however,  is  not  able  to  remove  interference  by 
higher  modes  which  are  superposed  on  the  fundmental  mode  in  the  correlation  functions. 

In  cases  of  modal  interference  or  multipathing,  modal  isolation  becomes  necessary. 
Time-variable  filtering  is  optimized  by  choosing  different  types  of  window  for  spectra  of 
different  degrees  of  sharpness,  and  by  selecting  the  proper  lags  for  the  window  used  in 
the  window-closing  procedure.  We  find  that  isolation  by  time-variable  filtering  is  better 
than  by  the  phase-matched  filtering  at  short  distances,  since  the  isolation  by  phase- 
matched  filtering  totally  depends  on  the  separation  between  the  signal  and  the  interfer¬ 
ing  noise  on  the  PAF.  After  isolation,  frequency-domain  Wiener  filtering  is  used  to  com¬ 
pute  the  interstation  Green’s  function  from  the  isolated  seismograms.  Even  though  the 
spectral  ratio  technique  will  work  almost  as  well  as  frequency-domain  Wiener  filtering 
when  there  is  little  noise  left  after  isolation,  the  frequency-domain  Wiener  filtering  can 
be  used  to  guarantee  the  removal  of  the  random  noise.  Group  and  phase  velocities,  and 


attenuation  coefficients  calculated  from  the  Green’s  function  by  the  proposed  technique 
are  very  stable  and  accurate.  In  addition,  the  MFT  is  found  to  have  large  spectral 
biases  which  result  in  errors  in  qr  estimates,  especially  for  the  higher  modes.  The  pro¬ 
posed  technique  in  this  paper  which  combines  modal  isolation  and  frequency-domain 
Wiener  filtering  works  well  for  higher  modes  as  well  as  for  the  fundamental  mode. 
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Synthetic  seismograms  and  spectra  of  Rayleigh  wave  generated  from  the 
model  listed  in  Table  1  at  distances  of  1000  km  and  2000  km  without  noise. 

Group  and  phase  velocities  and  attenuation  coefficients  by  different  tech¬ 
niques. 

The  seismograms  and  spectra  with  30%  Gaussian  random  noise  added  to  the 
seismogram  in  Figure  la.. 


l-' 


Windowing  on  correlation  functions  involved  in  time-domain  and  frequency- 
domain  Wiener  filtering.  Both  types  of  filtering  use  as  window’s  lag  for  the 
cross-correlation  function.  Frequency-domain  Wiener  filtering,  however,  uses 
a  window  of  smaller  lag,  ,  on  the  autocorrelation  function. 


Comparison  of  Green’s  functions  obtained  by  different  techniques. 


Comparison  of  the  group  and  phase  velocities  and  attenuation  coefficients  by 
different  techniques. 

Seismograms  and  spectra  composed  of  fundamental  and  first  higher  modes 
and  random  noise. 

Comparison  of  Green’s  functions  obtained  by  the  spectral  ratio  technique  and 
frequency-domain  Wiener  filtering. 


Cross-correlation  and  autocorrelation  functions  calculated  from  seismograms 
in  Figure  3a. 

The  seismogram,  amplitude  and  phase  spectra  at  2000  km  in  Figure  3a  before 
and  after  being  isolated  by  time-variable  filtering,  b  and  a  indicate  before 
and  and  after  the  isolation,  respectively. 


Comparison  of  group  and  phase  velocities  and  attenuation  coefficients 
obtained  by  frequency-domain  Wiener  filtering  after  isolation  by  time- 
variable  filtering  at  each  station  versus  theoretical  values. 

Comparison  of  Attenuation  coefficients  obtained  by  taking  the  spectral  ratio 
of  the  amplitudes  from  MFT  versus  theoretical  values. 


Comparison  of  spectral  amplitudes  and  attenuation  coefficients  of  the  higher 
mode  from  MFT  (circle)  and  from  time-variable  filtering  (dotted  solid  line) 
with  raw  spectra  (solid  line). 


The  pseudo- autocorrelation  function  in  phase-matclied  filtering,  and  the 
amplitude  spectrum  before  and  after  isolation  by  phase-matched  filtering  See 
Figure  4a  for  the  meaning  of  the  symbols  a  and  b. 

Synthetic  seismograms  and  spectra  at  2000  km  and  4000  km  composed  of 
fundamental  and  first  higher  modes  and  random  noise. 

Attenuation  coefficients  obtained  by  the  spectral  ratio  technique  without  and 
with  isolations  by  phase-matched  and  time-variable  filtering. 

The  vertical  component  seismograms  at  NDI  and  KOD  before  and  after  the 
isolation  by  time-variable  filtering.  See  Figure  4a  for  the  meaning  of  symbols 
a  and  b. 

Comparison  of  group  and  phase  velocities  and  attenuation  coefficients  by  the 
proposed  technique  and  the  spectral  ratio  technique. 
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Recently  developed  data  processing  methods  (Hwang  and  Mitchell, 
this  report)  and  newly  acquired  surface  wave  amplitude  data  have  allowed 
us  to  obtain  accurate  and  reliable  attenuation  coefficient  values  at 
periods  between  30  and  90  seconds  in  stable  continental  regions.  Previ¬ 
ous  studies,  such  as  Herrmann  and  Mitchell  (1975)  have  been  limited  by 


large  uncertainties  in  attenuation  coefficients  at  longer  periods. 
Because  of  those  uncertainties,  linear  inversions  for  Qp”^  produced 
poorly  resolved  models  with  very  low  values  (or  high  values)  in  the 
lower  crust  and  upper  mantle. 


The  new  attenuation  coefficient  values  extend  to  periods  as  great 
as  90  seconds  where  such  values  are  higher  than  those  previously 
obtained.  Figure  1  shows  an  example  for  eastern  South  America.  Note 
the  relatively  small  standard  deviations  compared  to  earlier  studies. 


Inversion  of  these  data  leads  to  a  model  with  relatively  low  values 
of  Qp  In  the  lower  crust  and  upper  mantle.  Figure  2  shows  the  model 
which  results  from  inverting  the  data  of  Figure  1 .  The  model  Is  charac¬ 
terized  by  high  Q  values  in  the  upper  crust  and  a  relatively  rapid  tran¬ 


sition  to  low  Q  values  at  mid-crustal  depths.  Note  that  this  transition 


occurs  at  depths  which  roughly  coincide  with  the  termination  of  Intra- 
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Figure  1.  Rayleigh  wave  attenuation  coefficients  obtained  by  Wiener  filtering 

and  modal  isolation  of  Interstatlon  data  across  eastern  South  Merlca 

Figure  2.  model  of  the  crust  and  i^per  mantle  beneath  eastern  South  America 

p 

obtained  by  inverting  the  attenuation  coefficient  data  of  Figure  1. 
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plate  earthquakes  In  stable  regions  and  may  represent  a  transition  from 
a  brittle  to  a  more  ductile  region  of  the  crust  (Brace  and  Kohlstedt, 
1980). 


Similar  results  have  been  obtained  for  the  Indian  shield  and  for 
recent  Inversions  of  new  data  In  the  eastern  United  States.  A  more  com¬ 
plete  presentation  of  the  Q  models  for  these  three  regions,  as  well  as 
tectonic  Implications  of  these  new  models,  will  be  presented  In  the  next 
technical  report. 
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